libint2-dev 2.3.0~beta3-2 / usr / include / libint2 / basis.h

/*
 *  This file is a part of Libint.
 *  Copyright (C) 2004-2014 Edward F. Valeev
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Library General Public License, version 2,
 *  as published by the Free Software Foundation.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU Library General Public License
 *  along with this program.  If not, see http://www.gnu.org/licenses/.
 *
 */

#ifndef _libint2_src_lib_libint_basis_h_
#define _libint2_src_lib_libint_basis_h_

#include <libint2/util/cxxstd.h>
#if LIBINT2_CPLUSPLUS_STD < 2011
# error "libint2/basis.h requires C++11 support"
#endif

#include <iostream>
#include <fstream>
#include <locale>
#include <vector>
#include <stdexcept>

#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include <libint2.h>
#include <libint2/shell.h>
#include <libint2/atom.h>

namespace libint2 {

  /// BasisSet is a slightly decorated \c std::vector of \c libint2::Shell objects.
  class BasisSet : public std::vector<libint2::Shell> {
    public:
      BasisSet() : name_(""), nbf_(-1), max_nprim_(0), max_l_(-1) {}
      BasisSet(const BasisSet&) = default;
      BasisSet(BasisSet&& other) :
        std::vector<libint2::Shell>(std::move(other)),
        name_(std::move(other.name_)),
        nbf_(other.nbf_),
        max_nprim_(other.max_nprim_),
        max_l_(other.max_l_),
        shell2bf_(std::move(other.shell2bf_))
      {
      }
      ~BasisSet() = default;
      BasisSet& operator=(const BasisSet&) = default;

      /// Construct from basis set name and vector of atoms

      /**
       * @param name basis set name
       * @param atoms \c std::vector of Atom objects
       * \note All atoms receive the same basis set
       * \note \c name will be "canonicalized" using BasisSet::canonicalize(name) to
       *       produce the file name where the basis will be sought. This file needs to contain
       *       the basis set definition in Gaussian94 format (see \c lib/basis directory for examples).
       *       The file will be searched in the following directories (in decreasing order of preference)
       *       <ol>
       *         <li> specified by LIBINT_DATA_PATH environmental variable, if defined </li>
       *         <li> specified by DATADIR macro variable, if defined </li>
       *         <li> specified by SRCDATADIR macro variable, if defined </li>
       *         <li> hardwired to directory \c /usr/share/libint/2.3.0/basis </li>
       *       </ol>
       */
      BasisSet(std::string name,
               const std::vector<Atom>& atoms) : name_(std::move(name)) {

        // read in the library file contents
        std::string basis_lib_path = data_path();
        auto canonical_name = canonicalize_name(name_);
        // some basis sets use cartesian d shells by convention, the convention is taken from Gaussian09
        auto force_cartesian_d = gaussian_cartesian_d_convention(canonical_name);
        
        // parse the name into components
        std::vector<std::string> basis_component_names = decompose_name_into_components(canonical_name);
        
        // ref_shells[component_idx][Z] => vector of Shells
        std::vector<std::vector<std::vector<libint2::Shell>>> component_basis_sets;
        component_basis_sets.reserve(basis_component_names.size());
        
        // read in ALL basis set components
        for(const auto& basis_component_name: basis_component_names) {
          auto file_dot_g94 = basis_lib_path + "/" + basis_component_name + ".g94";
        
          // use same cartesian_d convention for all components!
          component_basis_sets.emplace_back(read_g94_basis_library(file_dot_g94, force_cartesian_d));
        }
        
        // for each atom find the corresponding basis components
        for(auto a=0; a<atoms.size(); ++a) {

          auto Z = atoms[a].atomic_number;
          
          // add each component in order
          for(auto comp_idx=0; comp_idx!=component_basis_sets.size(); ++comp_idx) {
            const auto& component_basis_set = component_basis_sets[comp_idx];
            if (component_basis_set[Z].empty())
              throw std::string("did not find the basis for this Z in ") +
                    basis_lib_path + "/" + basis_component_names[comp_idx] + ".g94";

            for(auto s: component_basis_set[Z]) {
              this->push_back(std::move(s));
              this->back().move({{atoms[a].x, atoms[a].y, atoms[a].z}});
            } // shell loop
          } // basis component loop
        } // atom loop

        init();
      }

      /// forces solid harmonics/Cartesian Gaussians
      /// @param solid if true, force all shells with L>1 to be solid harmonics, otherwise force all shells to Cartesian
      void set_pure(bool solid) {
        for(auto& s: *this) {
          s.contr[0].pure = solid;
        }
        init();
      }

      /// @return the number of basis functions in the basis; -1 if uninitialized
      long nbf() const {
        return nbf_;
      }
      /// @return the maximum number of primitives in a contracted Shell, i.e. maximum contraction length; 0 if uninitialized
      size_t max_nprim() const {
        return max_nprim_;
      }
      /// @return the maximum angular momentum of a contraction; -1 if uninitialized
      long max_l() const {
        return max_l_;
      }
      /// @return the map from shell index to index of the first basis function from this shell
      /// \note basis functions are ordered as shells, i.e. shell2bf[i] >= shell2bf[j] iff i >= j
      const std::vector<size_t>& shell2bf() const {
        return shell2bf_;
      }
      /// Computes the map from shells to the corresponding atoms in \c atoms. If no atom matches the origin of a shell, it is mapped to -1.
      /// @note coordinates must match \em exactly , i.e. shell2atom[k] == l iff atoms[l].x == *this[k].O[0] && atoms[l].y == *this[k].O[1] &&  atoms[l].z == *this[k].O[2]
      /// @return the map from shell index to the atom in the list \c atoms that coincides with its origin;
      std::vector<long> shell2atom(const std::vector<Atom>& atoms) const {
        std::vector<long> result;
        result.reserve(size());
        for(const auto& s: *this) {
          auto a = std::find_if(atoms.begin(), atoms.end(), [&s](const Atom& a){ return s.O[0] == a.x && s.O[1] == a.y && s.O[2] == a.z; } );
          result.push_back( a != atoms.end() ? a - atoms.begin() : -1);
        }
        return result;
      }

    private:
      std::string name_;
      long nbf_;
      size_t max_nprim_;
      int max_l_;
      std::vector<size_t> shell2bf_;

      void init() {
        nbf_ = nbf(*this);
        max_nprim_ = max_nprim(*this);
        max_l_ = max_l(*this);
        shell2bf_ = compute_shell2bf(*this);
      }

      struct canonicalizer {
          char operator()(char c) {
            char cc = ::tolower(c);
            switch (cc) {
              case '/': cc = 'I'; break;
            }
            return cc;
          }
      };

      static std::string canonicalize_name(const std::string& name) {
        auto result = name;
        std::transform(name.begin(), name.end(),
                       result.begin(), BasisSet::canonicalizer());
        return result;
      }

      // see http://www.gaussian.com/g_tech/g_ur/m_basis_sets.htm
      bool gaussian_cartesian_d_convention(const std::string& canonical_name) {
        // 3-21??g??, 4-31g??
        if (canonical_name.find("3-21")    == 0 ||
            canonical_name.find("4-31g")   == 0)
          return true;
        // 6-31??g?? but not 6-311 OR 6-31g()
        if (canonical_name.find("6-31") == 0 && canonical_name[4] != '1') {
          // to exclude 6-31??g() find the g, then check the next character
          auto g_pos = canonical_name.find('g');
          if (g_pos == std::string::npos) // wtf, I don't even know what this is, assume spherical d is OK
            return false;
          if (g_pos+1 == canonical_name.size()) // 6-31??g uses cartesian d
            return true;
          if (canonical_name[g_pos+1] == '*') // 6-31??g*? uses cartesian d
            return true;
        }
        return false;
      }
      
      /// decompose basis set name into components
      std::vector<std::string> decompose_name_into_components(std::string name) {
        std::vector<std::string> component_names;
        // aug-cc-pvxz* = cc-pvxz* + augmentation-... , except aug-cc-pvxz-cabs
        if ( (name.find("aug-cc-pv") == 0) && (name.find("cabs")==std::string::npos)  ) {
          std::string base_name = name.substr(4);
          component_names.push_back(base_name);
          component_names.push_back(std::string("augmentation-") + base_name);
        }
        else
          component_names.push_back(name);

        return component_names;
      }
      
      /** determines the path to the data directory, as follows:
       *       <ol>
       *         <li> specified by LIBINT_DATA_PATH environmental variable, if defined </li>
       *         <li> specified by DATADIR macro variable, if defined </li>
       *         <li> specified by SRCDATADIR macro variable, if defined </li>
       *         <li> hardwired to directory \c /usr/share/libint/2.3.0/basis </li>
       *       </ol>
       *  @throw std::runtime_error if the path is not valid, or cannot be determined
       *  @return valid path to the data directory
       */
      static std::string data_path() {
        std::string path;
        const char* data_path_env = getenv("LIBINT_DATA_PATH");
        if (data_path_env) {
          path = data_path_env;
        }
        else {
#if defined(DATADIR)
          path = std::string{DATADIR};
#elif defined(SRCDATADIR)
          path = std::string{SRCDATADIR};
#else
          path = std::string("/usr/share/libint/2.3.0") + std::string("/basis");
#endif
        }
        // validate path
        bool error = true;
        if (not path.empty()) {
          struct stat sb;
          error = (::stat(path.c_str(), &sb) == -1);
          error = error && not S_ISDIR(sb.st_mode);
        }
        if (error) {
          std::ostringstream oss; oss << "BasisSet::data_path(): path \"" << path << "\" is not valid";
          throw std::runtime_error(oss.str());
        }
        return path;
      }
      
      /// converts fortran scientific-notation floats that use d/D instead of e/E in \c str
      /// @param[in,out] str string in which chars 'd' and 'D' are replaced with 'e' and 'E',
      ///                respectively
      static void fortran_dfloats_to_efloats(std::string& str) {
        for(auto& ch: str) {
          if (ch == 'd') ch = 'e';
          if (ch == 'D') ch = 'E';
        }
      }

    public:

      /** reads in all basis sets from a Gaussian94-formatted basis set file (see https://bse.pnl.gov/bse/portal)
       *  @param[in] file_dot_g94 file name
       *  @param[in] force_cartesian_d force use of Cartesian d shells, if true
       *  @param[in] locale_name specifies the locale to use
       *  @throw std::runtime_error if the path is not valid, or cannot be determined
       *  @return vector of basis sets for each element
       *  @warning the included library basis sets should be parsed using POSIX locale
       */
      static std::vector<std::vector<libint2::Shell>> read_g94_basis_library(std::string file_dot_g94,
                                                                             bool force_cartesian_d = false,
                                                                             std::string locale_name = std::string("POSIX")) {

        std::locale locale(locale_name.c_str());  // TODO omit c_str() with up-to-date stdlib
        std::cout << "Will read basis set from " << file_dot_g94 << std::endl;
        std::ifstream is(file_dot_g94);
        is.imbue(locale);
        if (not is.good()) {
          std::ostringstream oss;
          oss << "BasisSet::read_g94_basis_library(): could not open \"" << file_dot_g94 << "\"" << std::endl;
          throw std::runtime_error(oss.str());
        }
        std::vector<std::vector<libint2::Shell>> ref_shells(118); // 118 = number of chemical elements

        std::string line, rest;

        // skip till first basis
        while(std::getline(is, line) && line != "****") {
        }

#define LIBINT2_LINE_TO_STRINGSTREAM(line) \
        std::istringstream iss(line); \
        iss.imbue(locale);
        
        size_t Z;
        auto nextbasis = true, nextshell = false;
        // read lines till end
        while(std::getline(is, line)) {
          // skipping empties and starting with '!' (the comment delimiter)
          if (line.empty() || line[0] == '!') continue;
          if (line == "****") {
            nextbasis = true;
            nextshell = false;
            continue;
          }
          if (nextbasis) {
            nextbasis = false;
            LIBINT2_LINE_TO_STRINGSTREAM(line);
            std::string elemsymbol;
            iss >> elemsymbol >> rest;

            bool found = false;
            using libint2::chemistry::element_info;
            for(const auto& e: element_info) {
              if (strcaseequal(e.symbol, elemsymbol)) {
                Z = e.Z;
                found = true;
                break;
              }
            }
            if (not found) {
              std::ostringstream oss;
              oss << "in file " << file_dot_g94
                  << " found G94 basis set for element symbol \""
                  << elemsymbol << "\", not found in Periodic Table.";
              throw std::runtime_error(oss.str());
            }

            nextshell = true;
            continue;
          }
          if (nextshell) {
            LIBINT2_LINE_TO_STRINGSTREAM(line);
            std::string amlabel;
            unsigned nprim;
            iss >> amlabel >> nprim >> rest;
            if (amlabel != "SP" && amlabel != "sp") {
              assert(amlabel.size() == 1);
              auto l = Shell::am_symbol_to_l(amlabel[0]);
              std::vector<double> exps;
              std::vector<double> coeffs;
              for(auto p = 0; p!=nprim; ++p) {
                while(std::getline(is, line) && (line.empty() || line[0] == '!')) {}
                fortran_dfloats_to_efloats(line);
                LIBINT2_LINE_TO_STRINGSTREAM(line);
                double e, c;
                iss >> e >> c;
                exps.emplace_back(e);
                coeffs.emplace_back(c);
              }
              auto pure = force_cartesian_d ? (l>2) : (l>1);
              ref_shells[Z].push_back(
                  libint2::Shell{
                    exps,
                    {
                      {l, pure, coeffs}
                    },
                    {{0.0, 0.0, 0.0}}
                  }
              );
            }
            else { // split the SP shells
              std::vector<double> exps;
              std::vector<double> coeffs_s, coeffs_p;
              for(auto p = 0; p!=nprim; ++p) {
                while(std::getline(is, line) && (line.empty() || line[0] == '!')) {}
                fortran_dfloats_to_efloats(line);
                LIBINT2_LINE_TO_STRINGSTREAM(line);
                double e, c1, c2;
                iss >> e >> c1 >> c2;
                exps.emplace_back(e);
                coeffs_s.emplace_back(c1);
                coeffs_p.emplace_back(c2);
              }
              ref_shells[Z].push_back(
                  libint2::Shell{exps,
                {
                 {0, false, coeffs_s}
                },
                {{0.0, 0.0, 0.0}}
              }
              );
              ref_shells[Z].push_back(
                  libint2::Shell{ exps,
                {
                 {1, false, coeffs_p}
                },
                {{0.0, 0.0, 0.0}}
              }
              );
            }
          }
        }
        
#undef LIBINT2_LINE_TO_STRINGSTREAM
        
        return ref_shells;
      }

      static size_t nbf(const std::vector<libint2::Shell>& shells) {
        size_t n = 0;
        for (const auto& shell: shells)
          n += shell.size();
        return n;
      }

      static size_t max_nprim(const std::vector<libint2::Shell>& shells) {
        size_t n = 0;
        for (auto shell: shells)
          n = std::max(shell.nprim(), n);
        return n;
      }

      static int max_l(const std::vector<libint2::Shell>& shells) {
        int l = 0;
        for (auto shell: shells)
          for (auto c: shell.contr)
            l = std::max(c.l, l);
        return l;
      }

      static std::vector<size_t> compute_shell2bf(const std::vector<libint2::Shell>& shells) {
        std::vector<size_t> result;
        result.reserve(shells.size());

        size_t n = 0;
        for (auto shell: shells) {
          result.push_back(n);
          n += shell.size();
        }

        return result;
      }

  }; // BasisSet

} // namespace libint2

#endif /* _libint2_src_lib_libint_basis_h_ */