LIBINT  2.6.0
atom.h
1 /*
2  * Copyright (C) 2004-2019 Edward F. Valeev
3  *
4  * This file is part of Libint.
5  *
6  * Libint is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * Libint is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef _libint2_src_lib_libint_atom_h_
22 #define _libint2_src_lib_libint_atom_h_
23 
24 #include <libint2/util/cxxstd.h>
25 #if LIBINT2_CPLUSPLUS_STD < 2011
26 # error "libint2/atom.h requires C++11 support"
27 #endif
28 
29 #include <array>
30 #include <iostream>
31 #include <sstream>
32 #include <tuple>
33 #include <utility>
34 #include <vector>
35 
36 #include <libint2/chemistry/elements.h>
37 
38 namespace libint2 {
39 
40  struct Atom {
41  int atomic_number;
42  double x, y, z;
43  };
44  inline bool operator==(const Atom& atom1, const Atom& atom2) {
45  return atom1.atomic_number == atom2.atomic_number && atom1.x == atom2.x &&
46  atom1.y == atom2.y && atom1.z == atom2.z;
47  }
48 
49  namespace constants {
51  struct codata_2014 {
52  static constexpr double bohr_to_angstrom = 0.52917721067;
53  static constexpr double angstrom_to_bohr = 1 / bohr_to_angstrom;
54  };
56  struct codata_2010 {
57  static constexpr double bohr_to_angstrom = 0.52917721092;
58  static constexpr double angstrom_to_bohr = 1 / bohr_to_angstrom;
59  };
60  } // namespace constants
61 
62 } // namespace libint2
63 
64 namespace {
65 
66 bool strcaseequal(const std::string &a, const std::string &b) {
67  return a.size() == b.size() &&
68  std::equal(a.begin(), a.end(), b.begin(), [](char a, char b) {
69  return ::tolower(a) == ::tolower(b);
70  });
71 }
72 
76 inline std::tuple<std::vector<libint2::Atom>,
78 __libint2_read_dotxyz(std::istream &is, const double bohr_to_angstrom,
79  const bool pbc = false) {
80  using libint2::Atom;
81  const std::string caller = std::string("libint2::read_dotxyz") + (pbc ? "_pbc" : "");
82 
83  // first line = # of atoms
84  size_t natom;
85  is >> natom;
86 
87  // read off the rest of first line and discard
88  std::string rest_of_line;
89  std::getline(is, rest_of_line);
90 
91  // second line = comment
92  std::string comment;
93  std::getline(is, comment);
94 
95  // rest of lines are atoms (and unit cell parameters, if pbc = true)
96  const auto nlines_expected = natom + (pbc ? 3 : 0);
97  std::vector<Atom> atoms(natom, Atom{0, 0.0, 0.0, 0.0});
98  std::array<std::array<double, 3>, 3> unit_cell({{{0.0, 0.0, 0.0}}});
99  bool found_abc[3] = {false, false, false};
100  for (size_t line = 0, atom_index = 0; line < nlines_expected; ++line) {
101  if (is.eof())
102  throw std::logic_error(caller + ": expected " +
103  std::to_string(nlines_expected) +
104  " sets of coordinates but only " +
105  std::to_string(line) + " received");
106 
107  // read line
108  std::string linestr;
109  std::getline(is, linestr);
110  std::istringstream iss(linestr);
111  // then parse ... this handles "extended" XYZ formats
112  std::string element_symbol;
113  double x, y, z;
114  iss >> element_symbol >> x >> y >> z;
115 
116  // .xyz files report Cartesian coordinates in angstroms; convert to bohr
117  const auto angstrom_to_bohr = 1 / bohr_to_angstrom;
118 
119  auto assign_atom = [angstrom_to_bohr](Atom &atom, int Z, double x, double y,
120  double z) {
121  atom.atomic_number = Z;
122  atom.x = x * angstrom_to_bohr;
123  atom.y = y * angstrom_to_bohr;
124  atom.z = z * angstrom_to_bohr;
125  };
126  auto assign_xyz = [angstrom_to_bohr](std::array<double, 3> &xyz, double x,
127  double y, double z) {
128  xyz[0] = x * angstrom_to_bohr;
129  xyz[1] = y * angstrom_to_bohr;
130  xyz[2] = z * angstrom_to_bohr;
131  };
132 
133  auto axis = -1;
134  // if pbc = true, look for unit cell params
135  if (pbc) {
136  if (strcaseequal("AA", element_symbol))
137  axis = 0;
138  if (strcaseequal("BB", element_symbol))
139  axis = 1;
140  if (strcaseequal("CC", element_symbol))
141  axis = 2;
142  if (axis != -1) {
143  if (found_abc[axis] == true)
144  throw std::logic_error(
145  caller + ": unit cell parameter along Cartesian axis " +
146  std::to_string(axis) + " appears more than once");
147  assign_xyz(unit_cell[axis], x, y, z);
148  found_abc[axis] = true;
149  }
150  }
151 
152  // .xyz files report element labels, hence convert to atomic numbers
153  if (axis == -1) {
154  int Z = -1;
155  for (const auto &e : libint2::chemistry::get_element_info()) {
156  if (strcaseequal(e.symbol, element_symbol)) {
157  Z = e.Z;
158  break;
159  }
160  }
161  if (Z == -1) {
162  std::ostringstream oss;
163  oss << caller << ": element symbol \"" << element_symbol
164  << "\" is not recognized" << std::endl;
165  throw std::logic_error(oss.str().c_str());
166  }
167 
168  if (pbc && atom_index == atoms.size()) { // if PBC, check for too many atoms
169  throw std::logic_error(caller + ": too many atoms");
170  }
171  assign_atom(atoms[atom_index++], Z, x, y, z);
172  }
173  }
174 
175  // make sure all 3 axes were specified
176  if (pbc) {
177  for(auto xyz=0; xyz!=3; ++xyz)
178  if (found_abc[xyz] == false) {
179  throw std::logic_error(caller +
180  ": unit cell parameter along Cartesian axis " +
181  std::to_string(xyz) + " not given");
182  }
183  }
184 
185  return std::make_tuple(atoms, unit_cell);
186 }
187 
188 } // anonymous namespace
189 
190 namespace libint2 {
191 
205 inline std::vector<Atom> read_dotxyz(
206  std::istream &is,
207  const double bohr_to_angstrom = constants::codata_2010::bohr_to_angstrom) {
208  std::vector<Atom> atoms;
209  std::tie(atoms, std::ignore) =
210  __libint2_read_dotxyz(is, bohr_to_angstrom, false);
211  return atoms;
212 }
213 
229 inline auto read_dotxyz_pbc(
230  std::istream &is,
231  const double bohr_to_angstrom = constants::codata_2010::bohr_to_angstrom)
232  -> decltype(__libint2_read_dotxyz(is, bohr_to_angstrom, true)) {
233  return __libint2_read_dotxyz(is, bohr_to_angstrom, true);
234 }
235 
237 std::vector<std::pair<double, std::array<double, 3>>> inline make_point_charges(
238  const std::vector<libint2::Atom> &atoms) {
239  std::vector<std::pair<double, std::array<double, 3>>> q(atoms.size());
240  for (const auto &atom : atoms) {
241  q.emplace_back(static_cast<double>(atom.atomic_number),
242  std::array<double, 3>{{atom.x, atom.y, atom.z}});
243  }
244  return q;
245 }
246 
247 } // namespace libint2
248 
249 #endif /* _libint2_src_lib_libint_atom_h_ */
std::vector< std::pair< double, std::array< double, 3 > > > make_point_charges(const std::vector< libint2::Atom > &atoms)
converts a vector of Atoms to a vector of point charges
Definition: atom.h:237
the 2010 CODATA reference set, available at DOI 10.1103/RevModPhys.84.1527
Definition: atom.h:56
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
Definition: atom.h:40
the 2014 CODATA reference set, available at DOI 10.1103/RevModPhys.88.035009
Definition: atom.h:51
auto read_dotxyz_pbc(std::istream &is, const double bohr_to_angstrom=constants::codata_2010::bohr_to_angstrom) -> decltype(__libint2_read_dotxyz(is, bohr_to_angstrom, true))
reads the list of atoms from a file in the PBC-extended XYZ format
Definition: atom.h:229
std::vector< Atom > read_dotxyz(std::istream &is, const double bohr_to_angstrom=constants::codata_2010::bohr_to_angstrom)
reads the list of atoms from a file in the standard XYZ format supported by most chemistry software (...
Definition: atom.h:205
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14