basisset.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2020 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
// VOTCA includes
#include <votca/tools/property.h>
 
// Local VOTCA includes
#include "votca/xtp/basisset.h"
#include <votca/tools/globals.h>
 
namespace votca {
namespace xtp {
 
L StringToEnum(const std::string& type) {
  assert(!type.empty() && "Shelltype must be non empty!");
  assert(type.size() == 1 && "Shelltype size must be one");
  const char t = type.back();
  return StringToEnum(t);
}
 
L StringToEnum(char type) {
  L l;
  if (type == 'S') {
    l = L::S;
  } else if (type == 'P') {
    l = L::P;
  } else if (type == 'D') {
    l = L::D;
  } else if (type == 'F') {
    l = L::F;
  } else if (type == 'G') {
    l = L::G;
  } else if (type == 'H') {
    l = L::H;
  } else if (type == 'I') {
    l = L::I;
  } else {
    throw std::runtime_error("FindLmax: Shelltype '" + std::string(1, type) +
                             "' not known");
  }
  return l;
}
 
std::string EnumToString(L l) {
  switch (l) {
    case L::S:
      return "S";
    case L::P:
      return "P";
    case L::D:
      return "D";
    case L::F:
      return "F";
    case L::G:
      return "G";
    case L::H:
      return "H";
    case L::I:
      return "I";
  }
  return "";
}
 
Index OffsetFuncShell(L l) {
  switch (l) {
    case L::S:
      return 0;
    case L::P:
      return 1;
    case L::D:
      return 4;
    case L::F:
      return 9;
    case L::G:
      return 16;
    case L::H:
      return 25;
    case L::I:
      return 36;
  }
  return -1;
}
 
Index NumFuncShell(L l) { return 2 * Index(l) + 1; }
 
Index NumFuncShell_cartesian(L l) {
  Index lindex = Index(l);
  return (lindex + 1) * (lindex + 2) / 2;
}
 
Index OffsetFuncShell_cartesian(L l) {
  switch (l) {
    case L::S:
      return 0;
    case L::P:
      return 1;
    case L::D:
      return 4;
    case L::F:
      return 10;
    case L::G:
      return 20;
    case L::H:
      return 35;
    case L::I:
      return 56;
  }
  return -1;
}
 
bool CheckShellType(const std::string& shelltype) {
  if (shelltype.empty()) {
    return false;
  }
  std::vector<char> allowed_shells = {'S', 'P', 'D', 'F', 'G', 'H', 'I'};
  std::vector<char>::iterator it =
      std::find(allowed_shells.begin(), allowed_shells.end(), shelltype[0]);
  if (it == allowed_shells.end()) {
    return false;
  } else {
    Index index = std::distance(allowed_shells.begin(), it);
    for (Index i = 1; i < Index(shelltype.size()); i++) {
      if (index + i > Index(allowed_shells.size()) ||
          shelltype[i] != allowed_shells[index + i]) {
        return false;
      }
    }
  }
 
  return true;
}
 
void BasisSet::Load(const std::string& name) {
 
  // if name contains .xml, assume a basisset .xml file is located in the
  // working directory
  std::string xmlFile;
  if (name.find(".xml") != std::string::npos) {
    xmlFile = name;
  } else {
    xmlFile = tools::GetVotcaShare() + "/xtp/basis_sets/" + name + ".xml";
  }
  tools::Property basis_property;
  basis_property.LoadFromXML(xmlFile);
  name_ = basis_property.get("basis").getAttribute<std::string>("name");
  std::vector<tools::Property*> elementProps =
      basis_property.Select("basis.element");
 
  for (tools::Property* elementProp : elementProps) {
    std::string elementName = elementProp->getAttribute<std::string>("name");
    Element& element = addElement(elementName);
    std::vector<tools::Property*> shellProps = elementProp->Select("shell");
    for (tools::Property* shellProp : shellProps) {
      std::string shellType = shellProp->getAttribute<std::string>("type");
      if (!CheckShellType(shellType)) {
        throw std::runtime_error("Shelltype: '" + shellType +
                                 "' is not a valid shelltype!");
      }
      for (char subtype : shellType) {
 
        double shellScale = shellProp->getAttribute<double>("scale");
 
        Shell& shell = element.addShell(StringToEnum(subtype), shellScale);
        std::vector<tools::Property*> constProps =
            shellProp->Select("constant");
        for (tools::Property* constProp : constProps) {
          double decay = constProp->getAttribute<double>("decay");
          std::vector<tools::Property*> contrProps =
              constProp->Select("contractions");
          double contraction = 0.0;
          for (tools::Property* contrProp : contrProps) {
            std::string contrType =
                contrProp->getAttribute<std::string>("type");
            if (contrType != std::string(1, subtype)) {
              continue;
            }
            contraction = contrProp->getAttribute<double>("factor");
          }
          shell.addGaussian(decay, contraction);
        }
      }
    }
  }
  return;
}
 
// adding an Element to a Basis Set
Element& BasisSet::addElement(std::string elementType) {
  auto e = elements_.insert({elementType, Element(elementType)});
  if (!e.second) {
    throw std::runtime_error("Inserting element into basisset failed!");
  }
  return e.first->second;
}
 
const Element& BasisSet::getElement(std::string element_type) const {
  std::map<std::string, Element>::const_iterator itm =
      elements_.find(element_type);
  if (itm == elements_.end()) {
    throw std::runtime_error("Basis set " + name_ +
                             " does not have element of type " + element_type);
  }
  return itm->second;
}
 
std::ostream& operator<<(std::ostream& out, const Shell& shell) {
 
  out << "Type:" << EnumToString(shell.getL()) << " Scale:" << shell.getScale()
      << " Func: " << shell.getnumofFunc() << "\n";
  for (const auto& gaussian : shell.gaussians_) {
    out << " Gaussian Decay: " << gaussian.decay();
    out << " Contraction: " << gaussian.contraction();
    out << "\n";
  }
  return out;
}
 
std::ostream& operator<<(std::ostream& out, const Element& element) {
  out << "Element:" << element.getType() << "\n";
  for (const auto& shell : element) {
    out << shell;
  }
  return out;
}
 
std::ostream& operator<<(std::ostream& out, const BasisSet& basis) {
  out << "BasisSet:" << basis.name_ << "\n";
  for (const auto& element : basis) {
    out << element.second;
  }
  out << std::flush;
  return out;
}
 
GaussianPrimitive& Shell::addGaussian(double decay, double contraction) {
  gaussians_.push_back(GaussianPrimitive(decay, contraction));
  return gaussians_.back();
}
 
}  // namespace xtp
}  // namespace votca