psml.hpp

Go to the documentation of this file.

#ifndef PSEUDO_PSML_HPP
#define PSEUDO_PSML_HPP
 
/*
 Copyright (C) 2018 Xavier Andrade
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 3 of the License, or
 (at your option) any later version.
 
 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 Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/
 
#include <cassert>
#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
 
#include "anygrid.hpp"
#include "base.hpp"
#include "element.hpp"
#include <rapidxml.hpp>
 
namespace pseudopotential {
 
class psml : public pseudopotential::anygrid {
 
public:
  psml(const std::string &filename, bool uniform_grid = false)
      : pseudopotential::anygrid(uniform_grid), file_(filename.c_str()),
        buffer_((std::istreambuf_iterator<char>(file_)),
                std::istreambuf_iterator<char>()) {
 
    filename_ = filename;
 
    buffer_.push_back('\0');
    doc_.parse<0>(&buffer_[0]);
 
    root_node_ = doc_.first_node("psml");
 
    spec_node_ = root_node_->first_node("pseudo-atom-spec");
    // some files do not have "pseudo-atom-spec" but "header"
    if (!spec_node_)
      spec_node_ = root_node_->first_node("header");
    assert(spec_node_);
 
    // now check the type
    bool has_local_potential = root_node_->first_node("local-potential");
    bool has_nl_projectors = select_set_node("nonlocal-projectors");
    bool has_semilocal_potentials = select_set_node("semilocal-potentials");
    bool has_pseudo_wavefunctions = select_set_node("pseudo-wave-functions");
    if (has_nl_projectors && has_local_potential) {
      type_ = pseudopotential::type::KLEINMAN_BYLANDER;
    } else if (has_semilocal_potentials && has_pseudo_wavefunctions) {
      type_ = pseudopotential::type::SEMILOCAL;
    } else {
      throw status::UNSUPPORTED_TYPE;
    }
 
    {
      // read lmax
      std::string tag1, tag2;
      if (type_ == pseudopotential::type::KLEINMAN_BYLANDER) {
        tag1 = "nonlocal-projectors";
        tag2 = "proj";
      } else if (type_ == pseudopotential::type::SEMILOCAL) {
        tag1 = "semilocal-potentials";
        tag2 = "slps";
      } else {
        assert(false);
      }
 
      lmax_ = -1;
      rapidxml::xml_node<> *node = select_set_node(tag1.c_str());
      assert(node);
      node = node->first_node(tag2.c_str());
      while (node) {
        int read_l = letter_to_l(node->first_attribute("l")->value());
        lmax_ = std::max(lmax_, read_l);
        node = node->next_sibling(tag2.c_str());
      }
      assert(lmax_ >= 0);
      assert(lmax_ < 9);
    }
 
    // read grid
    {
      rapidxml::xml_node<> *node = root_node_->first_node("grid");
 
      assert(node);
 
      int size = value<int>(node->first_attribute("npts"));
      grid_.resize(size);
      std::istringstream stst(node->first_node("grid-data")->value());
      for (int ii = 0; ii < size; ii++) {
        stst >> grid_[ii];
        if (ii > 0)
          assert(grid_[ii] > grid_[ii - 1]);
      }
 
      assert(fabs(grid_[0]) <= 1e-10);
 
      mesh_size_ = 0;
      for (double rr = 0.0; rr <= grid_[grid_.size() - 1]; rr += mesh_spacing())
        mesh_size_++;
 
      grid_weights_.resize(grid_.size());
 
      // the integration weights are not available, we approximate them
      grid_weights_[0] = 0.5 * (grid_[1] - grid_[0]);
      for (unsigned ii = 1; ii < grid_.size() - 1; ii++)
        grid_weights_[ii] = grid_[ii + 1] - grid_[ii - 1];
      grid_weights_[grid_.size() - 1] =
          0.5 * (grid_[grid_.size() - 1] - grid_[grid_.size() - 2]);
    }
  }
 
  pseudopotential::format format() const {
    return pseudopotential::format::PSML;
  }
 
  int size() const { return buffer_.size(); };
 
  std::string description() const { return ""; }
 
  std::string symbol() const {
    return element::trim(spec_node_->first_attribute("atomic-label")->value());
  }
 
  int atomic_number() const {
    return value<int>(spec_node_->first_attribute("atomic-number"));
  }
 
  double mass() const {
    element el(symbol());
    return el.mass();
  }
 
  double valence_charge() const {
    return value<double>(spec_node_->first_node("valence-configuration")
                             ->first_attribute("total-valence-charge"));
  }
 
  int llocal() const { return -1; }
 
  pseudopotential::exchange exchange() const {
    // PSML uses libxc ids, so we just need to read the value
    rapidxml::xml_node<> *node = spec_node_->first_node("exchange-correlation")
                                     ->first_node("libxc-info")
                                     ->first_node("functional");
    while (node) {
      if (value<std::string>(node->first_attribute("type")) == "exchange") {
        return pseudopotential::exchange(
            value<int>(node->first_attribute("id")));
      }
      node = node->next_sibling("functional");
    }
    return pseudopotential::exchange::UNKNOWN;
  }
 
  pseudopotential::correlation correlation() const {
    // PSML uses libxc ids, so we just need to read the value
    rapidxml::xml_node<> *node = spec_node_->first_node("exchange-correlation")
                                     ->first_node("libxc-info")
                                     ->first_node("functional");
    while (node) {
      if (value<std::string>(node->first_attribute("type")) == "correlation") {
        return pseudopotential::correlation(
            value<int>(node->first_attribute("id")));
      }
      node = node->next_sibling("functional");
    }
    return pseudopotential::correlation::UNKNOWN;
  }
 
  int nchannels() const {
    if (type_ == pseudopotential::type::SEMILOCAL)
      return 1;
    if (is_fully_relativistic()) {
      int nc = 0;
      for (int l = 0; l <= lmax_; l++)
        nc = std::max(nc, nprojectors_per_l(l));
      return nc;
    }
 
    int nc = 0;
    rapidxml::xml_node<> *node = nonlocal_projectors_node();
    assert(node);
    node = node->first_node("proj");
    while (node) {
      int read_ic = value<int>(node->first_attribute("seq")) - 1;
      nc = std::max(nc, read_ic + 1);
      node = node->next_sibling("proj");
    }
    return nc;
  }
 
  void local_potential(std::vector<double> &val) const {
    read_function(root_node_->first_node("local-potential"), val, true);
  }
 
  int nprojectors() const {
    rapidxml::xml_node<> *node =
        nonlocal_projectors_node()->first_node("proj");
    int count = 0;
    while (node) {
      count++;
      node = node->next_sibling("proj");
    }
    return count;
  }
 
  int nprojectors_per_l(int l) const {
    if (type_ == pseudopotential::type::SEMILOCAL)
      return 1;
 
    rapidxml::xml_node<> *node =
        nonlocal_projectors_node()->first_node("proj");
    int count = 0;
    while (node) {
      int read_l = letter_to_l(node->first_attribute("l")->value());
      if (read_l == l) count++; 
      node = node->next_sibling("proj");
    } 
    return count;
  }
 
 
  void projector(int l, int ic, std::vector<double> &val) const {
    rapidxml::xml_node<> *node = projector_node(l, ic);
    read_function(node, val);
  }
 
  double d_ij(int l, int ic, int jc) const {
    if (ic != jc)
      return 0.0;
 
    rapidxml::xml_node<> *node = projector_node(l, ic);
    if (!node)
      return 0.0;
    return value<double>(node->first_attribute("ekb"));
  }
 
  bool has_radial_function(int l) const {
    return wavefunction_node_for_l(l) != nullptr;
  }
 
  void radial_function(int l, std::vector<double> &val) const {
    rapidxml::xml_node<> *node = wavefunction_node_for_l(l);
    read_function(node, val);
  }
 
  void radial_potential(int l, std::vector<double> &val) const {
    rapidxml::xml_node<> *node =
        semilocal_potentials_node()->first_node("slps");
    while (node) {
      int read_l = letter_to_l(node->first_attribute("l")->value());
      if (l == read_l)
        break;
      node = node->next_sibling("slps");
    }
    read_function(node, val);
  }
 
  bool has_nlcc() const {
    rapidxml::xml_node<> *node = root_node_->first_node("pseudocore-charge");
    return node;
  }
 
  void nlcc_density(std::vector<double> &val) const {
    read_function(root_node_->first_node("pseudocore-charge"), val);
    for (unsigned ii = 0; ii < val.size(); ii++)
      val[ii] /= 4.0 * M_PI;
  }
 
  bool has_density() { return root_node_->first_node("valence-charge"); }
 
  void density(std::vector<double> &val) const {
    read_function(root_node_->first_node("valence-charge"), val);
    for (unsigned ii = 0; ii < val.size(); ii++)
      val[ii] /= 4.0 * M_PI;
  }
 
  bool has_total_angular_momentum() const {
    return spec_node_->first_attribute("relativity")->value() ==
           std::string("dirac");
  }
 
  int projector_2j(int l, int ic) const {
    rapidxml::xml_node<> *node = projector_node(l, ic);
    assert(node);
    double read_j = value<double>(node->first_attribute("j"));
    return lrint(2.0 * read_j);
  }
 
  int nwavefunctions() const {
    rapidxml::xml_node<> *wf_set = pseudo_wavefunctions_node();
    if (!wf_set)
      return 0;
 
    int count = 0;
    for (rapidxml::xml_node<> *shell = first_valence_shell();
         shell; shell = shell->next_sibling("shell")) {
      rapidxml::xml_node<> *node = wf_set->first_node("pswf");
      while (node) {
        if (wavefunction_matches_shell(node, shell))
          count++;
        node = node->next_sibling("pswf");
      }
    }
 
    if (count > 0)
      return count;
 
    count = 0;
    rapidxml::xml_node<> *node = wf_set->first_node("pswf");
    while (node) {
      count++;
      node = node->next_sibling("pswf");
    }
    return count;
  }
 
  void wavefunction(int index, int &n, int &l, double &occ,
                    std::vector<double> &val) const {
    rapidxml::xml_node<> *node = wavefunction_node(index);
    assert(node);
 
    n = value<int>(node->first_attribute("n"));
    l = letter_to_l(node->first_attribute("l")->value());
    occ = wavefunction_occupation(node);
 
    read_function(node, val);
  }
 
  int wavefunction_2j(int ii) const {
    rapidxml::xml_node<> *node = wavefunction_node(ii - 1);
    assert(node);
 
    auto *j_attr = node->first_attribute("j");
    if (!j_attr)
      return 0;
 
    return lrint(2.0 * value<double>(j_attr));
  }
 
private:
  bool is_fully_relativistic() const {
    auto *attr = spec_node_->first_attribute("relativity");
    return attr && std::string(attr->value()) == "dirac";
  }
 
  rapidxml::xml_node<> *select_set_node(const char *name) const {
    rapidxml::xml_node<> *node = root_node_->first_node(name);
    rapidxml::xml_node<> *fallback = nullptr;
    const std::string preferred_set =
        is_fully_relativistic() ? "lj" : "scalar_relativistic";
 
    while (node) {
      auto *set_attr = node->first_attribute("set");
      if (!set_attr)
        return node;
 
      std::string set_name = set_attr->value();
      if (set_name == preferred_set)
        return node;
      if (set_name == "scalar_relativistic" && fallback == nullptr)
        fallback = node;
      if (fallback == nullptr)
        fallback = node;
 
      node = node->next_sibling(name);
    }
 
    return fallback;
  }
 
  rapidxml::xml_node<> *nonlocal_projectors_node() const {
    return select_set_node("nonlocal-projectors");
  }
 
  rapidxml::xml_node<> *semilocal_potentials_node() const {
    return select_set_node("semilocal-potentials");
  }
 
  rapidxml::xml_node<> *pseudo_wavefunctions_node() const {
    return select_set_node("pseudo-wave-functions");
  }
 
  rapidxml::xml_node<> *projector_node(int l, int ic) const {
    assert(l >= 0);
    rapidxml::xml_node<> *node = nonlocal_projectors_node()->first_node("proj");
    int current_ic = 0;
    int desired_2j = 0;
    int desired_seq = 0;
 
    if (is_fully_relativistic() && l > 0) {
      desired_seq = ic / 2 + 1;
      desired_2j = (ic % 2 == 0) ? 2 * l + 1 : 2 * l - 1;
    }
 
    while (node) {
      int read_l = letter_to_l(node->first_attribute("l")->value());
      if (l == read_l) {
        if (is_fully_relativistic()) {
          if (l == 0) {
            if (current_ic == ic)
              return node;
            current_ic++;
          } else {
            int read_ic = value<int>(node->first_attribute("seq"));
            int read_2j = lrint(2.0 * value<double>(node->first_attribute("j")));
            if (read_ic == desired_seq && read_2j == desired_2j)
              return node;
          }
        } else {
          int read_ic = value<int>(node->first_attribute("seq")) - 1;
          if (ic == read_ic)
            return node;
        }
      }
      node = node->next_sibling("proj");
    }
 
    return nullptr;
  }
 
  rapidxml::xml_node<> *first_valence_shell() const {
    rapidxml::xml_node<> *node = spec_node_->first_node("valence-configuration");
    if (!node)
      return nullptr;
    return node->first_node("shell");
  }
 
  bool wavefunction_matches_shell(rapidxml::xml_node<> *wavefunction,
                                  rapidxml::xml_node<> *shell) const {
    assert(wavefunction);
    assert(shell);
 
    return value<int>(wavefunction->first_attribute("n")) ==
               value<int>(shell->first_attribute("n")) &&
           std::string(wavefunction->first_attribute("l")->value()) ==
               shell->first_attribute("l")->value();
  }
 
  double wavefunction_occupation(rapidxml::xml_node<> *wavefunction) const {
    assert(wavefunction);
 
    auto *occupation_attr = wavefunction->first_attribute("occupation");
    if (occupation_attr)
      return value<double>(occupation_attr);
 
    rapidxml::xml_node<> *shell = first_valence_shell();
    while (shell) {
      if (wavefunction_matches_shell(wavefunction, shell)) {
        double occupation = value<double>(shell->first_attribute("occupation"));
        auto *j_attr = wavefunction->first_attribute("j");
        if (!j_attr)
          return occupation;
 
        double j = value<double>(j_attr);
        int l = letter_to_l(wavefunction->first_attribute("l")->value());
        return occupation * (2.0 * j + 1.0) / (2.0 * (2.0 * l + 1.0));
      }
      shell = shell->next_sibling("shell");
    }
 
    return 0.0;
  }
 
  rapidxml::xml_node<> *wavefunction_node_for_l(int l) const {
    rapidxml::xml_node<> *wf_set = pseudo_wavefunctions_node();
    if (!wf_set)
      return nullptr;
 
    rapidxml::xml_node<> *node = wf_set->first_node("pswf");
    while (node) {
      int read_l = letter_to_l(node->first_attribute("l")->value());
      if (l == read_l) {
        rapidxml::xml_node<> *shell = first_valence_shell();
        while (shell) {
          if (wavefunction_matches_shell(node, shell))
            return node;
          shell = shell->next_sibling("shell");
        }
        return node;
      }
      node = node->next_sibling("pswf");
    }
 
    return nullptr;
  }
 
  rapidxml::xml_node<> *wavefunction_node(int index) const {
    rapidxml::xml_node<> *wf_set = pseudo_wavefunctions_node();
    assert(wf_set);
 
    int current = 0;
    for (rapidxml::xml_node<> *shell = first_valence_shell();
         shell; shell = shell->next_sibling("shell")) {
      rapidxml::xml_node<> *best_high_j = nullptr;
      rapidxml::xml_node<> *best_low_j = nullptr;
 
      for (rapidxml::xml_node<> *node = wf_set->first_node("pswf"); node;
           node = node->next_sibling("pswf")) {
        if (!wavefunction_matches_shell(node, shell))
          continue;
 
        auto *j_attr = node->first_attribute("j");
        if (!j_attr) {
          if (current == index)
            return node;
          current++;
          continue;
        }
 
        double j = value<double>(j_attr);
        int shell_l = letter_to_l(shell->first_attribute("l")->value());
        if (j > shell_l)
          best_high_j = node;
        else
          best_low_j = node;
      }
 
      if (best_high_j) {
        if (current == index)
          return best_high_j;
        current++;
      }
      if (best_low_j) {
        if (current == index)
          return best_low_j;
        current++;
      }
    }
 
    for (rapidxml::xml_node<> *node = wf_set->first_node("pswf"); node;
         node = node->next_sibling("pswf")) {
      if (current == index)
        return node;
      current++;
    }
 
    return nullptr;
  }
 
  void read_function(rapidxml::xml_node<> *base_node, std::vector<double> &val,
                     bool potential_padding = false) const {
    assert(base_node);
    rapidxml::xml_node<> *node =
        base_node->first_node("radfunc")->first_node("data");
    assert(node);
    int size = grid_.size();
    if (node->first_attribute("npts"))
      size = value<int>(node->first_attribute("npts"));
    val.resize(grid_.size());
    std::istringstream stst(node->value());
    for (int ii = 0; ii < size; ii++)
      stst >> val[ii];
 
    if (potential_padding) {
      for (unsigned ii = size; ii < grid_.size(); ii++)
        val[ii] = -valence_charge() / grid_[ii];
    } else {
      for (unsigned ii = size; ii < grid_.size(); ii++)
        val[ii] = 0.0;
    }
 
    interpolate(val);
  }
 
  // for some stupid reason psml uses letters instead of numbers for angular
  // momentum
  static int letter_to_l(const std::string &letter) {
    if (letter == "s")
      return 0;
    if (letter == "p")
      return 1;
    if (letter == "d")
      return 2;
    if (letter == "f")
      return 3;
    if (letter == "g")
      return 4;
    if (letter == "h")
      return 5;
    return -1;
  }
 
  std::ifstream file_;
  std::vector<char> buffer_;
  rapidxml::xml_document<> doc_;
  rapidxml::xml_node<> *root_node_;
  rapidxml::xml_node<> *spec_node_;
};
 
} // namespace pseudopotential
 
#endif