LIBINT  2.6.0
molden.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 INCLUDE_LIBINT2_LCAO_MOLDEN_H_
22 #define INCLUDE_LIBINT2_LCAO_MOLDEN_H_
23 
24 #include <iomanip>
25 #include <ostream>
26 #include <string>
27 #include <vector>
28 
29 #include <libint2/atom.h>
30 #include <libint2/basis.h>
31 #include <libint2/cgshell_ordering.h>
32 #include <libint2/chemistry/elements.h>
33 #include <libint2/shell.h>
34 #include <libint2/shgshell_ordering.h>
35 
36 #pragma GCC diagnostic push
37 #pragma GCC system_header
38 #include <Eigen/Core>
39 #pragma GCC diagnostic pop
40 
41 namespace libint2 {
42 namespace molden {
43 
46 class Export {
47  public:
73  template <typename Coeffs, typename Occs, typename Energies = Eigen::VectorXd>
74  Export(const std::vector<Atom>& atoms, const std::vector<Shell>& basis,
75  const Coeffs& coefficients, const Occs& occupancies,
76  const Energies& energies = Energies(),
77  const std::vector<std::string>& symmetry_labels =
78  std::vector<std::string>(),
79  const std::vector<bool>& spincases = std::vector<bool>(),
80  double coefficient_epsilon = 5e-11)
81  : atoms_(atoms),
82  basis_(validate(basis)),
83  coefs_(coefficients),
84  occupancies_(occupancies),
85  energies_(energies),
86  labels_(symmetry_labels),
87  spins_(spincases),
88  coefficient_epsilon_(coefficient_epsilon) {
89  initialize_bf_map();
90  }
91 
93  void write_prologue(std::ostream& os) const {
94  os << "[Molden Format]" << std::endl;
95  }
96 
98  void write_atoms(std::ostream& os) const {
99  os << "[Atoms] AU" << std::endl;
100 
101  os.fill(' ');
102  os << std::fixed << std::setprecision(8);
103  auto iatom = 0;
104  for (const auto& atom : atoms_) {
105  auto Z = atom.atomic_number;
106  os << std::setw(4) << libint2::chemistry::get_element_info().at(Z - 1).symbol
107  << std::setw(6) << (iatom + 1) << std::setw(6) << Z << std::setw(14)
108  << atom.x << std::setw(14) << atom.y << std::setw(14) << atom.z
109  << std::endl;
110  ++iatom;
111  }
112  }
113 
116  void write_basis(std::ostream& os) const {
117  bool f_found = false;
118  os << "[GTO]" << std::endl;
119  for (size_t iatom = 0; iatom < atoms_.size(); ++iatom) {
120  os << std::setw(4) << (iatom + 1) << std::setw(4) << 0 << std::endl;
121  for (auto ish : atom2shell_[iatom]) {
122  const Shell& sh = basis_.at(ish);
123  if (sh.contr.size() == 1) {
124  const auto& contr = sh.contr[0];
125  const auto l = contr.l;
126  assert(l <= 4); // only up to g functions are supported
127  if (l == 3)
128  f_found = true;
129  const auto nprim = contr.coeff.size();
130  os << std::setw(4) << Shell::am_symbol(contr.l) << std::setw(6)
131  << nprim << std::setw(6) << "1.00" << std::endl;
132  for (int iprim = 0; iprim < nprim; ++iprim) {
133  os << std::scientific << std::uppercase << std::setprecision(10);
134  os << std::setw(20) << sh.alpha[iprim] << std::setw(20)
135  << sh.coeff_normalized(0, iprim) << std::endl;
136  } // end loop over primitives
137  } // end if ncontraction == 1
138  else {
139  assert(false); // Not implemented
140  }
141  } // end loop over shells on center
142  // format calls for a blank line here
143  os << std::endl;
144  } // end loop over centers
145 
146  // write solid harmonic/cartesian tags
147  {
148  // Molden default is cartesians throughout
149  // dfg_is_cart_ is set to true even if there are no shells of a given type
150  if (dfg_is_cart_[0]) { // cartesian d
151  if (!dfg_is_cart_[1]) // solid harmonic f
152  os << "[7F]" << std::endl;
153  } else { // solid harmonic d
154  if (!dfg_is_cart_[1]) { // solid harmonic f
155  os << "[5D7F]" << std::endl;
156  } else if (f_found) { // cartesian f
157  os << "[5D10F]" << std::endl;
158  } else { // no f functions
159  os << "[5D]" << std::endl;
160  }
161  }
162  if (!dfg_is_cart_[2]) // solid harmonic g
163  os << "[9G]" << std::endl;
164  }
165  }
166 
168  void write_lcao(std::ostream& os) const {
169  os << "[MO]" << std::endl;
170  for (int imo = 0; imo < coefs_.cols(); ++imo) {
171  os << std::fixed << std::setprecision(10);
172  os << std::setw(8) << "Sym= " << (labels_.empty() ? "A" : labels_.at(imo))
173  << std::endl
174  << std::setw(8) << "Ene= " << std::setw(16)
175  << (energies_.rows() == 0 ? 0.0 : energies_(imo)) << std::endl
176  << std::setw(8) << "Spin= "
177  << (spins_.empty() ? "Alpha" : (spins_.at(imo) ? "Alpha" : "Beta"))
178  << std::endl
179  << std::setw(8) << "Occup= " << occupancies_(imo) << std::endl;
180  os << std::scientific << std::uppercase << std::setprecision(10);
181  for (int iao = 0; iao < coefs_.rows(); ++iao) {
182  const auto C_ao_mo = coefs_(ao_map_[iao], imo);
183  if (std::abs(C_ao_mo) >= coefficient_epsilon_) {
184  os << std::setw(6) << (iao + 1) << " " << std::setw(16)
185  << C_ao_mo << std::endl;
186  }
187  } // end loop over AOs
188  } // end loop over MOs
189  }
190 
192  void write(std::ostream& os) const {
193  write_prologue(os);
194  write_atoms(os);
195  write_basis(os);
196  write_lcao(os);
197  }
198 
200  void write(const std::string& filename) const {
201  std::ofstream os(filename);
202  write_prologue(os);
203  write_atoms(os);
204  write_basis(os);
205  write_lcao(os);
206  }
207 
208  private:
209  const std::vector<Atom>& atoms_;
210  const std::vector<Shell>& basis_;
211  Eigen::MatrixXd coefs_;
212  Eigen::VectorXd occupancies_;
213  Eigen::VectorXd energies_;
214  std::vector<std::string> labels_;
215  std::vector<bool> spins_;
216  double coefficient_epsilon_;
217  mutable bool
218  dfg_is_cart_[3]; // whether {d, f, g} shells are cartesian (true) or
219  // solid harmonics (false)
220  std::vector<std::vector<long>>
221  atom2shell_; // maps atom -> shell indices in basis_
222  std::vector<long>
223  ao_map_; // maps from the AOs ordered according to Molden
224  // (atoms->shells, bf in shells ordered in the Molden order)
225  // to the AOs ordered according to basis_
226 
229  const std::vector<Shell>& validate(const std::vector<Shell>& shells) const {
230  bool dfg_found[] = {false, false, false};
231  for(int i=0; i!=sizeof(dfg_is_cart_)/sizeof(bool); ++i)
232  dfg_is_cart_[i] = true;
233  for (const auto& shell : shells) {
234  for (const auto& contr : shell.contr) {
235  if (contr.l > 4)
236  throw std::logic_error(
237  "molden::Export cannot handle shells with l > 4");
238 
239  switch (contr.l) {
240  case 1:
241  if (contr.pure)
242  throw std::logic_error(
243  "molden::Export cannot handle solid harmonics p shells");
244  break;
245  case 2:
246  case 3:
247  case 4: {
248  if (!dfg_found[contr.l - 2]) {
249  dfg_is_cart_[contr.l - 2] = !contr.pure;
250  dfg_found[contr.l - 2] = true;
251  }
252  if (!contr.pure ^ dfg_is_cart_[contr.l - 2])
253  throw std::logic_error(
254  "molden::Export only supports all-Cartesian or "
255  "all-solid-harmonics d/f/g shells");
256  }
257 
258  default: {} // l = 0 is fine
259  }
260  }
261  }
262  return shells;
263  }
264 
265  void initialize_bf_map() {
266  atom2shell_ = BasisSet::atom2shell(atoms_, basis_);
267 
268  const auto nao = BasisSet::nbf(basis_);
269  ao_map_.resize(nao);
270  assert(nao == coefs_.rows());
271  const auto shell2ao = BasisSet::compute_shell2bf(basis_);
272  long ao_molden = 0;
273  for (size_t iatom = 0; iatom < atoms_.size(); ++iatom) {
274  for (auto ish : atom2shell_[iatom]) {
275  auto ao = shell2ao[ish]; // refers to order assumed by coefs
276  const auto& shell = basis_[ish];
277  const auto ncontr = shell.contr.size();
278  for (int c = 0; c != ncontr; ++c) {
279  const auto l = shell.contr[c].l;
280  const auto pure = shell.contr[c].pure;
281  if (pure) {
282  int m;
283  FOR_SOLIDHARM_MOLDEN(l, m)
284  const auto ao_in_shell = INT_SOLIDHARMINDEX(l, m);
285  ao_map_[ao_molden] = ao + ao_in_shell;
286  ++ao_molden;
287  END_FOR_SOLIDHARM_MOLDEN
288  ao += 2 * l + 1;
289  } else {
290  int i, j, k;
291  FOR_CART_MOLDEN(i, j, k, l)
292  const auto ao_in_shell = INT_CARTINDEX(l, i, j);
293  ao_map_[ao_molden] = ao + ao_in_shell;
294  ++ao_molden;
295  END_FOR_CART_MOLDEN
296  ao += INT_NCART(l);
297  }
298  } // contraction loop
299  }
300  }
301  }
302 
303 }; // Export
304 
307 class PBCExport: public Export{
308  public:
333  template <typename Coeffs, typename Occs, typename Energies = Eigen::VectorXd>
334  PBCExport(const std::vector<Atom>& atoms,
335  const Eigen::VectorXd& cell_axes,
336  const std::vector<Shell>& basis,
337  const Coeffs& coefficients,
338  const Occs& occupancies,
339  int space_group,
340  const Energies& energies = Energies(),
341  const std::vector<std::string>& symmetry_labels =
342  std::vector<std::string>(),
343  const std::vector<bool>& spincases = std::vector<bool>())
344  : Export(atoms, basis, coefficients, occupancies, energies, symmetry_labels, spincases),
345  cell_axes_(cell_axes),
346  space_group_(space_group)
347  {
348  //initialize_bf_map();
349  }
350 
352  void write_space_group(std::ostream& os) const {
353  os << "[SpaceGroup] (Number)" << std::endl;
354  os << space_group_ << std::endl;
355  }
356 
358  void write_operators(std::ostream& os) const {
359  os << "[Operators]" << std::endl;
360  os << "x, y, z" << std::endl;
361  }
362 
364  void write_cell_axes(std::ostream& os) const {
365 
366  os << "[CellAxes] (Ang)" << std::endl;
367  os << std::setw(4) << cell_axes_[0] << std::setw(12) << 0.0 << std::setw(12) << 0.0 << std::endl;
368  os << std::setw(4) << 0.0 << std::setw(12) << cell_axes_[1] << std::setw(12) << 0.0 << std::endl;
369  os << std::setw(4) << 0.0 << std::setw(12) << 0.0 << std::setw(12) << cell_axes_[2] << std::endl;
370  }
371 
372  void write(const std::string& filename) const {
373  std::ofstream os(filename);
374  write_prologue(os);
375  write_space_group(os);
376  write_operators(os);
377  write_cell_axes(os);
378  write_atoms(os);
379  write_basis(os);
380  write_lcao(os);
381  }
382 
383 private:
384  Eigen::VectorXd cell_axes_;
385  int space_group_;
386 
387 }; // PBCExport
388 
389 } // namespace molden
390 } // namespace libint2
391 
392 #endif // INCLUDE_LIBINT2_LCAO_MOLDEN_H_
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:81
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
svector< Contraction > contr
contractions
Definition: shell.h:105
void write(std::ostream &os) const
writes "prologue", atoms, basis, and LCAOs to ostream os
Definition: molden.h:192
void write_lcao(std::ostream &os) const
writes the "[MO]" section to ostream os
Definition: molden.h:168
Export(const std::vector< Atom > &atoms, const std::vector< Shell > &basis, const Coeffs &coefficients, const Occs &occupancies, const Energies &energies=Energies(), const std::vector< std::string > &symmetry_labels=std::vector< std::string >(), const std::vector< bool > &spincases=std::vector< bool >(), double coefficient_epsilon=5e-11)
Definition: molden.h:74
void write_operators(std::ostream &os) const
writes the "[Operators]" section to ostream os
Definition: molden.h:358
void write_cell_axes(std::ostream &os) const
writes the "[CellAxes]" section to ostream os
Definition: molden.h:364
void write_prologue(std::ostream &os) const
writes "[Molden Format]" to ostream os
Definition: molden.h:93
real_t coeff_normalized(size_t c, size_t p) const
Definition: shell.h:225
Exports LCAO coefficients in Molden format.
Definition: molden.h:46
PBCExport(const std::vector< Atom > &atoms, const Eigen::VectorXd &cell_axes, const std::vector< Shell > &basis, const Coeffs &coefficients, const Occs &occupancies, int space_group, const Energies &energies=Energies(), const std::vector< std::string > &symmetry_labels=std::vector< std::string >(), const std::vector< bool > &spincases=std::vector< bool >())
Definition: molden.h:334
svector< real_t > alpha
exponents
Definition: shell.h:104
Extension of the Molden exporter to support JMOL extensions for crystal orbitals (see https://sourcef...
Definition: molden.h:307
void write(const std::string &filename) const
same as write(ostream), but creates new file named filename
Definition: molden.h:200
void write_atoms(std::ostream &os) const
writes the "[Atoms]" section to ostream os
Definition: molden.h:98
void write_space_group(std::ostream &os) const
writes the "[SpaceGroup]" section to ostream os
Definition: molden.h:352
void write_basis(std::ostream &os) const
writes the "[GTO]" section, as well as optional Cartesian/solid harmonics keywords,...
Definition: molden.h:116