molecule.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.
 
  Copyright (C) 2007,2010 Oak Ridge National Laboratory
 
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 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 General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
  For more information please contact:
 
  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367
 
  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680
*/
 
#ifndef MADNESS_CHEM_MOLECULE_H__INCLUDED
#define MADNESS_CHEM_MOLECULE_H__INCLUDED
 
/// \file moldft/molecule.h
/// \brief Declaration of molecule related classes and functions
 
#include <ctype.h>
#include <madness/chem/atomutil.h>
#include <madness/chem/corepotential.h>
#include <madness/misc/misc.h>
#include <madness/tensor/tensor.h>
#include <madness/world/vector.h>
 
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
 
#include "madness/mra/QCCalculationParametersBase.h"
#include "madness/mra/commandlineparser.h"
 
namespace madness {
 
class World;
 
class Atom {
 public:
  double x, y, z, q;           ///< Coordinates and charge in atomic units
  unsigned int atomic_number;  ///< Atomic number
  double mass;                 ///< Mass
  bool pseudo_atom;            ///< Indicates if this atom uses a pseudopotential
 
  explicit Atom(double x, double y, double z, double q, unsigned int atomic_number, bool pseudo_atom)
      : x(x), y(y), z(z), q(q), atomic_number(atomic_number), pseudo_atom(pseudo_atom) {
    mass = get_atomic_data(atomic_number).mass;
 
    if (mass == -1.0) MADNESS_EXCEPTION("faulty element in Atom", 1);
 
    // unstable elements are indicated by negative masses, the mass
    // is taken from the longest-living element
    if (mass < 0.0) mass *= -1.0;
  }
 
  explicit Atom(double x, double y, double z, double q, unsigned int atomic_number) : x(x), y(y), z(z), q(q), atomic_number(atomic_number) {
    mass = get_atomic_data(atomic_number).mass;
 
    if (mass == -1.0) MADNESS_EXCEPTION("faulty element in Atom", 1);
 
    // unstable elements are indicated by negative masses, the mass
    // is taken from the longest-living element
    if (mass < 0.0) mass *= -1.0;
 
    pseudo_atom = false;
  }
 
  Atom(const Atom &a) : x(a.x), y(a.y), z(a.z), q(a.q), atomic_number(a.atomic_number), mass(a.mass), pseudo_atom(a.pseudo_atom) {}
 
  /// Default construct makes a zero charge ghost atom at origin
  Atom() : x(0), y(0), z(0), q(0), atomic_number(0), mass(0.0), pseudo_atom(false) {}
 
  bool operator==(const Atom &other) const {
    double thresh=1.e-10;
    auto displacement=Vector<double, 3>({x, y, z})-(Vector<double, 3>({other.x, other.y, other.z}));
    double err = displacement.normf();
    return ((err<thresh) && q == other.q &&
            atomic_number == other.atomic_number && mass == other.mass &&
            pseudo_atom == other.pseudo_atom);
  }
 
  int get_atomic_number() const { return atomic_number; }
 
  madness::Vector<double, 3> get_coords() const { return madness::Vector<double, 3>{x, y, z}; }
 
  /// return the mass in atomic units (electron mass = 1 a.u.)
  double get_mass_in_au() const { return constants::atomic_mass_in_au * mass; }
 
  template <typename Archive>
  void serialize(Archive &ar) {
    ar & x & y & z & q & atomic_number & mass & pseudo_atom;
  }
  hashT hash() const {
    hashT h = hash_value(x);
    hash_combine(h, y);
    hash_combine(h, z);
    hash_combine(h, q);
    hash_combine(h, atomic_number);
    hash_combine(h, mass);
    hash_combine(h, pseudo_atom);
    return h;
  }
};
 
std::ostream &operator<<(std::ostream &s, const Atom &atom);
 
class Molecule {
 public:
  // Needed for ParameterManager
  static constexpr char const *tag = "molecule";
  [[nodiscard]] json to_json_if_precedence(std::string const &precedence) const {
    json mol_schema = to_json();
    mol_schema["parameters"] = parameters.to_json_if_precedence(precedence);
    insert_symbols_and_geometry(mol_schema);
    return mol_schema;
  }
 
 
  struct GeometryParameters : public QCCalculationParametersBase {
    GeometryParameters(const GeometryParameters &other) = default;
 
    GeometryParameters(World &world, const commandlineparser &parser) : GeometryParameters() {
      try {
        set_global_convenience_options(parser);
        // The deck has to carry a geometry group -- unless the geometry already
        // came in through `--geometry=`/`--molecule=`, in which case there is
        // nothing to read and a missing group is not an error. This used to be
        // unconditionally fatal and was papered over by
        // read_input_and_commandline_options swallowing the exception.
        throw_if_datagroup_not_found = not is_user_defined("source_name");
        // read the data group the input file actually uses, then apply
        // `--geometry=...`/`--molecule=...` overrides for the other spelling
        // as well (a no-op if that key is absent)
        const std::string tag = input_tag(parser);
        read_input_and_commandline_options(world, parser, tag);
        read_commandline_options(world, parser, tag == "geometry" ? "molecule" : "geometry");
        set_derived_values(parser);
 
      } catch (std::exception &e) {
        print(tag, "end");
        throw;
        //                MADNESS_EXCEPTION("faulty geometry input",1);
      }
    }
 
    GeometryParameters() {
      ignore_unknown_keys = true;
      ignore_unknown_keys_silently = true;
      throw_if_datagroup_not_found = true;
 
      //            initialize<std::vector<std::string>>("source",{"inputfile"},"where
      //            to get the coordinates from: ({inputfile}, {library,xxx},
      //            {xyz,xxx.xyz})");
      initialize<std::string>("source_type", "inputfile", "where to get the coordinates from", {"inputfile", "xyz", "library"});
      initialize<std::string>("source_name", "TBD", "name of the geometry from the library or the input file");
      initialize<double>("eprec", 1.e-4, "smoothing for the nuclear potential");
      initialize<std::string>("units", "atomic", "coordinate units", {"atomic", "angstrom", "bohr", "au"});
      initialize<std::vector<double>>("field", {0.0, 0.0, 0.0}, "external electric field");
      initialize<bool>("no_orient", false,
                       "if true the coordinates will not be reoriented and/or symmetrized");
      initialize<double>("symtol", -1.e-2,
                         "distance threshold");
 
      initialize<std::string>("core_type", "none", "core potential type", {"none", "mcp"});
      initialize<bool>("psp_calc", false, "pseudopotential calculation for all atoms");
      initialize<bool>("pure_ae", true, "pure all electron calculation with no pseudo-atoms");
    }
 
    std::string get_tag() const override {
      return std::string("molecule");
    }
 
    /// source_name is a path (xyz file / input file) -- never case-fold it
    std::set<std::string> case_sensitive_keys() const override {
      return {"source_name"};
    }
 
    void set_global_convenience_options(const commandlineparser &parser) {
      // `--geometry=<name>` is the documented spelling (see the help text of
      // moldft/nemo/cc2/madqc); `--molecule=<name>` is accepted as an alias.
      // value_raw, not value: this may be a path, and paths are case-sensitive.
      for (const char *key : {"geometry", "molecule"}) {
        if (parser.key_exists(key)) {
          set_user_defined_value("source_name", parser.value_raw(key));
        }
      }
    }
 
    /// name of the data group carrying the geometry in the input file
 
    /// The group is spelled `geometry` (canonical -- the moldft/nemo/cc2 decks)
    /// or `molecule` (alias -- the madqc decks and the structure library).
    /// Returns whichever the input file contains, preferring the canonical name
    /// and falling back to it when the file has neither (or does not exist), so
    /// that the missing-data-group diagnostic names the documented spelling.
    static std::string input_tag(const commandlineparser &parser);
 
    void set_derived_values(const commandlineparser &parser) {
      // check if we use an xyz file, the structure library or the input file
      set_derived_value("source_name",
                        parser.value("input"));  // will not override user input
      std::string src_type = derive_source_type_from_name(source_name(), parser);
      set_derived_value("source_type", src_type);
      if (parser.key_exists("no_orient") and parser.value("no_orient") == "true") set_derived_value("no_orient", true);
 
      // check for ambiguities in the derived source type
      if (not is_user_defined("source_type")) {
        std::ifstream f(source_name().c_str());
        bool found_geometry_file = f.good();
        //                bool
        //                found_geometry_file=std::filesystem::exists(source_name());
 
        bool geometry_found_in_library = true;
        try {  // check for existence of file and structure in the library
          std::ifstream f;
          position_stream_in_library(f, source_name());
        } catch (...) {
          geometry_found_in_library = false;
        }
 
        if (found_geometry_file and geometry_found_in_library) {
          madness::print("\n\n");
          madness::print(
              "geometry specification ambiguous: found geometry in "
              "the structure library and in a file\n");
          madness::print("  ", get_structure_library_path());
          madness::print("  ", source_name());
          madness::print(
              "\nPlease specify the location of your geometry input "
              "by one of the two lines:\n");
          madness::print("  source_type xyz");
          madness::print("  source_type library\n\n");
          MADNESS_EXCEPTION("faulty input\n\n", 1);
        }
      }
 
      //            std::vector<std::string> src=source();
      //
      //            // some convenience for the user
      //
      //            // if source is the input file provide the name of the input
      //            file if (src.size()==1 and src[0]=="inputfile")
      //                set_derived_value("source",std::vector<std::string>({src[0],parser.value("input")}));
      //            // if source is not "inputfile" or "library" assume an xyz
      //            file if (src.size()==1 and src[0]!="inputfile") {
      //                std::size_t found=src[0].find("xyz");
      //                if (found==src[0].size()-3) { // check input file ends
      //                with xyz
      //                    set_user_defined_value("source",
      //                    std::vector<std::string>({"xyz", src[0]}));
      //                } else {
      //                    throw std::runtime_error("error in deriving geometry
      //                    parameters");
      //                }
      //            }
 
      // Resolve a relative xyz path to an absolute one *now*, while the process
      // is still in the directory the user invoked it from: madqc chdir's into a
      // per-task run directory before the task re-reads its geometry, and it
      // propagates source_name through the deck it writes there, so a relative
      // path would then resolve against the task directory and not be found.
      // Only "xyz" is rewritten -- for "inputfile" the name must keep comparing
      // equal to parser.value("input") in derive_source_type_from_name(), and
      // "library" names are not paths.
      if (source_type() == "xyz") {
        const std::string abs = absolute_path(source_name());
        if (abs != source_name()) set_user_defined_value("source_name", abs);
      }
 
      if (source_type() == "xyz") set_derived_value("units", std::string("angstrom"));
      if (units() == "bohr" or units() == "au") set_derived_value("units", std::string("atomic"));
    }
 
    /// make a relative path absolute, relative to the current working directory
 
    /// Returns \p path unchanged if it is already absolute or cannot be resolved
    /// here -- the reader then reports the name as given, which is the more
    /// useful diagnostic.
    static std::string absolute_path(const std::string &path);
 
    std::string source_type() const { return get<std::string>("source_type"); }
    std::string source_name() const { return get<std::string>("source_name"); }
    std::vector<double> field() const { return get<std::vector<double>>("field"); }
    double eprec() const { return get<double>("eprec"); }
    std::string units() const { return get<std::string>("units"); }
    std::string core_type() const { return get<std::string>("core_type"); }
    bool psp_calc() const { return get<bool>("psp_calc"); }
    bool pure_ae() const { return get<bool>("pure_ae"); }
    bool no_orient() const { return get<bool>("no_orient"); }
    double symtol() const { return get<double>("symtol"); }
 
    static std::string derive_source_type_from_name(const std::string name, const commandlineparser &parser) {
      if (name == parser.value("input")) return "inputfile";
      std::size_t pos = name.find(".xyz");
      if (pos != std::string::npos) return "xyz";
      return "library";
    }
  };
 
 private:
  // If you add more fields don't forget to serialize them
  std::vector<Atom> atoms;
  std::vector<double> rcut;  // Reciprocal of the smoothing radius
  CorePotentialManager core_pot;
  madness::Tensor<double> field;
 
  /// The molecular point group is automatically assigned in the
  /// identify_pointgroup function
  std::string pointgroup_ = "c1";
 
 public:
  GeometryParameters parameters;
 
  static void print_parameters();
 
  std::string get_pointgroup() const { return pointgroup_; }
 
 private:
  void swapaxes(int ix, int iy);
 
  template <typename opT>
  bool test_for_op(opT op, const double symtol) const;
 
  template <typename opT>
  void symmetrize_for_op(opT op, const double symtol);
 
  template <typename opT>
  int find_symmetry_equivalent_atom(int iatom, opT op, const double symtol) const;
 
  bool test_for_c2(double xaxis, double yaxis, double zaxis, const double symtol) const;
 
  bool test_for_sigma(double xaxis, double yaxis, double zaxis, const double symtol) const;
 
  bool test_for_inverse(const double symtol) const;
 
  /// Apply to (x,y,z) a C2 rotation about an axis thru the origin and
  /// (xaxis,yaxis,zaxis)
  struct apply_c2 {
    double xaxis, yaxis, zaxis;
    apply_c2(double xaxis, double yaxis, double zaxis) : xaxis(xaxis), yaxis(yaxis), zaxis(zaxis) {}
    void operator()(double &x, double &y, double &z) const {
      double raxissq = xaxis * xaxis + yaxis * yaxis + zaxis * zaxis;
      double dx = x * xaxis * xaxis / raxissq;
      double dy = y * yaxis * yaxis / raxissq;
      double dz = z * zaxis * zaxis / raxissq;
      x = 2.0 * dx - x;
      y = 2.0 * dy - y;
      z = 2.0 * dz - z;
    }
  };
 
  /// Apply to (x,y,z) a reflection through a plane containing the origin with
  /// normal (xaxis,yaxis,zaxis)
  struct apply_sigma {
    double xaxis, yaxis, zaxis;
    apply_sigma(double xaxis, double yaxis, double zaxis) : xaxis(xaxis), yaxis(yaxis), zaxis(zaxis) {}
    void operator()(double &x, double &y, double &z) const {
      double raxissq = xaxis * xaxis + yaxis * yaxis + zaxis * zaxis;
      double dx = x * xaxis * xaxis / raxissq;
      double dy = y * yaxis * yaxis / raxissq;
      double dz = z * zaxis * zaxis / raxissq;
 
      x = x - 2.0 * dx;
      y = y - 2.0 * dy;
      z = z - 2.0 * dz;
    }
  };
 
  struct apply_inverse {
    double xaxis, yaxis, zaxis;
    apply_inverse(double xaxis, double yaxis, double zaxis) : xaxis(xaxis), yaxis(yaxis), zaxis(zaxis) {}
    void operator()(double &x, double &y, double &z) const {
      x = -x;
      y = -y;
      z = -z;
    }
  };
 
 public:
  /// Makes a molecule with zero atoms
  Molecule() : atoms(), rcut(), core_pot(), field(3L){};
 
  /// makes a molecule from a list of atoms
  Molecule(std::vector<Atom> atoms, double eprec, CorePotentialManager core_pot = {}, madness::Tensor<double> field = madness::Tensor<double>(3L));
 
  /// makes a molecule using contents of \p parser
  Molecule(World &world, const commandlineparser &parser);
 
  void get_structure();
 
  void read_structure_from_library(const std::string &name);
 
  static std::istream &position_stream_in_library(std::ifstream &f, const std::string &name);
 
  static std::string get_structure_library_path();
 
  /// print out a Gaussian cubefile header
 
  /// @param[in] offset  the offset to be subtracted from the coordinates
  std::vector<std::string> cubefile_header(const Vector<double, 3> offset = Vector<double, 3>(0.0)) const;
 
  // initializes Molecule using the contents of file \c filename
  void read_file(const std::string &filename);
 
  // initializes Molecule using the contents of stream \c f
  void read(std::istream &f);
 
  // initializes Molecule using the contents of file \c filenam assuming an xyz
  // file
  void read_xyz(const std::string filename);
 
  void read_core_file(const std::string &filename);
 
  std::string guess_file() const { return core_pot.guess_file(); };
 
  unsigned int n_core_orb_all() const;
 
  unsigned int n_core_orb(unsigned int atn) const {
    if (core_pot.is_defined(atn))
      return core_pot.n_core_orb_base(atn);
    else
      return 0;
  };
 
  unsigned int get_core_l(unsigned int atn, unsigned int c) const { return core_pot.get_core_l(atn, c); }
 
  double get_core_bc(unsigned int atn, unsigned int c) const { return core_pot.get_core_bc(atn, c); }
 
  double core_eval(int atom, unsigned int core, int m, double x, double y, double z) const;
 
  double core_derivative(int atom, int axis, unsigned int core, int m, double x, double y, double z) const;
 
  bool is_potential_defined(unsigned int atn) const { return core_pot.is_defined(atn); };
 
  bool is_potential_defined_atom(int i) const { return core_pot.is_defined(atoms[i].atomic_number); };
 
  void add_atom(double x, double y, double z, double q, int atn);
 
  void add_atom(double x, double y, double z, double q, int atn, bool psat);
 
  size_t natom() const { return atoms.size(); };
 
  void set_atom_charge(unsigned int i, double zeff);
 
  unsigned int get_atom_charge(unsigned int i) const;
 
  unsigned int get_atomic_number(unsigned int i) const;
 
  void set_pseudo_atom(unsigned int i, bool psat);
 
  bool get_pseudo_atom(unsigned int i) const;
 
  void set_atom_coords(unsigned int i, double x, double y, double z);
 
  madness::Tensor<double> get_all_coords() const;
 
  std::vector<madness::Vector<double, 3>> get_all_coords_vec() const;
 
  std::vector<double> atomic_radii;
 
  void set_all_coords(const madness::Tensor<double> &newcoords);
 
  void update_rcut_with_eprec(double value);
 
  void set_rcut(double value);
 
  std::vector<double> get_rcut() const { return rcut; }
 
  void set_core_eprec(double value) { core_pot.set_eprec(value); }
 
  void set_core_rcut(double value) { core_pot.set_rcut(value); }
 
  double get_eprec() const { return parameters.eprec(); }
 
  double bounding_cube() const;
 
  const Atom &get_atom(unsigned int i) const;
 
  const std::vector<Atom> &get_atoms() const { return atoms; }
 
  /// print all parameters and the molecular geometry
  void print(std::ostream& os=std::cout, const bool defined_only=false) const;
 
  /// print user-defined parameters and the molecular geometry
  void print_defined_only(std::ostream& os=std::cout) const {
      this->print(os,true);
  };
 
  double inter_atomic_distance(unsigned int i, unsigned int j) const;
 
  double nuclear_repulsion_energy() const;
 
  double nuclear_repulsion_derivative(size_t iatom, int axis) const;
 
  /// compute the nuclear-nuclear contribution to the second derivatives
 
  /// @param[in]  iatom   the i-th atom (row of the hessian)
  /// @param[in]  jatom   the j-th atom (column of the hessian)
  /// @param[in]  iaxis   the xyz axis of the i-th atom
  /// @param[in]  jaxis   the xyz axis of the j-th atom
  /// return the (3*iatom + iaxis, 3*jatom + jaxis) matix element of the hessian
  double nuclear_repulsion_second_derivative(int iatom, int jatom, int iaxis, int jaxis) const;
 
  /// return the hessian matrix of the second derivatives d^2/dxdy V
 
  /// no factor 0.5 included
  Tensor<double> nuclear_repulsion_hessian() const;
 
  /// compute the dipole moment of the nuclei
 
  ///  @param[in] axis the axis (x, y, z)
  double nuclear_dipole(int axis) const;
 
  /// compute the derivative of the nuclear dipole wrt a nuclear displacement
 
  /// @param[in]  atom    the atom which will be displaced
  /// @param[in]  axis    the axis where the atom will be displaced
  /// @return     a vector which all 3 components of the dipole derivative
  Tensor<double> nuclear_dipole_derivative(const int atom, const int axis) const;
 
  /// evaluate the nuclear charge density at point `{x,y,z}` using the default
  /// MADNESS nuclear model. See smoothed_density() for the description
  /// of the default nuclear model.
  /// \param x,y,z the point at which the nuclear charge density is evaluated
  /// \param rscale setting `rscale>1` will make a nucleus larger by a factor of
  /// \p rscale (in other words, `rcut` is multiplied by the inverse of by this)
  /// \return the nuclear charge density at point `{x,y,z}`
  /// \sa smoothed_density()
  double nuclear_charge_density(double x, double y, double z, double rscale = 1.) const;
 
  double smallest_length_scale() const;
 
  std::string symmetrize_and_identify_point_group(const double symtol);
 
  /// Moves the center of nuclear charge to the origin
  void center();
 
  /// rotates the molecule and the external field
 
  void fix_phase();
 
  /// @param[in]  D   the rotation matrix
  void rotate(const Tensor<double> &D);
 
  /// translate the molecule
  void translate(const Tensor<double> &translation);
 
  Tensor<double> center_of_mass() const;
 
  /// compute the mass-weighting matrix for the hessian
 
  /// use as
  /// mass_weighted_hessian=inner(massweights,inner(hessian,massweights));
  Tensor<double> massweights() const {
    Tensor<double> M(3 * natom(), 3 * natom());
    for (size_t i = 0; i < natom(); i++) {
      const double sqrtmass = 1.0 / sqrt(get_atom(i).get_mass_in_au());
      M(3 * i, 3 * i) = sqrtmass;
      M(3 * i + 1, 3 * i + 1) = sqrtmass;
      M(3 * i + 2, 3 * i + 2) = sqrtmass;
    }
    return M;
  }
 
  Tensor<double> moment_of_inertia() const;
 
  void orient(bool verbose = false);
 
  double total_nuclear_charge() const;
 
  /// nuclear attraction potential for the whole molecule
  double nuclear_attraction_potential(double x, double y, double z) const;
 
  /// nuclear attraction potential for a specific atom in the molecule
  double atomic_attraction_potential(int iatom, double x, double y, double z) const;
 
  double molecular_core_potential(double x, double y, double z) const;
 
  double core_potential_derivative(int atom, int axis, double x, double y, double z) const;
 
  double nuclear_attraction_potential_derivative(int atom, int axis, double x, double y, double z) const;
 
  double nuclear_attraction_potential_second_derivative(int atom, int iaxis, int jaxis, double x, double y, double z) const;
 
 bool operator==(const Molecule& other) const {
    if (atoms.size() != other.atoms.size() || rcut.size() != other.rcut.size() || pointgroup_ != other.pointgroup_) return false;
    for (size_t i = 0; i < atoms.size(); ++i) {
      if (not (atoms[i] == other.atoms[i])) return false;
    }
    for (size_t i = 0; i < rcut.size(); ++i) {
      if (rcut[i] != other.rcut[i]) return false;
    }
 
    return (field-other.field).normf()<1.e-13 && parameters == other.parameters;
  }
 
  template <typename Archive>
  void serialize(Archive &ar) {
    ar & atoms & rcut & core_pot & parameters & pointgroup_ & field;
  }
 
  hashT hash() const {
    hashT h = hash_range(atoms.begin(), atoms.end());
    hash_combine(h, hash_range(rcut.begin(), rcut.end()));
    hash_combine(h, pointgroup_);
    return h;
  }
  [[nodiscard]] json to_json() const;
  void insert_symbols_and_geometry(json &mol_json) const;
  void from_json(const json &mol_json);
};
 
}  // namespace madness
 
#endif