21 #ifndef INCLUDE_LIBINT2_LCAO_MOLDEN_H_ 22 #define INCLUDE_LIBINT2_LCAO_MOLDEN_H_ 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> 36 #pragma GCC diagnostic push 37 #pragma GCC system_header 39 #pragma GCC diagnostic pop 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)
82 basis_(validate(basis)),
84 occupancies_(occupancies),
86 labels_(symmetry_labels),
88 coefficient_epsilon_(coefficient_epsilon) {
94 os <<
"[Molden Format]" << std::endl;
99 os <<
"[Atoms] AU" << std::endl;
102 os << std::fixed << std::setprecision(8);
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
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;
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)
150 if (dfg_is_cart_[0]) {
151 if (!dfg_is_cart_[1])
152 os <<
"[7F]" << std::endl;
154 if (!dfg_is_cart_[1]) {
155 os <<
"[5D7F]" << std::endl;
156 }
else if (f_found) {
157 os <<
"[5D10F]" << std::endl;
159 os <<
"[5D]" << std::endl;
162 if (!dfg_is_cart_[2])
163 os <<
"[9G]" << std::endl;
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))
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"))
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;
192 void write(std::ostream& os)
const {
200 void write(
const std::string& filename)
const {
201 std::ofstream os(filename);
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_;
220 std::vector<std::vector<long>>
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) {
236 throw std::logic_error(
237 "molden::Export cannot handle shells with l > 4");
242 throw std::logic_error(
243 "molden::Export cannot handle solid harmonics p shells");
248 if (!dfg_found[contr.l - 2]) {
249 dfg_is_cart_[contr.l - 2] = !contr.pure;
250 dfg_found[contr.l - 2] =
true;
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");
265 void initialize_bf_map() {
266 atom2shell_ = BasisSet::atom2shell(atoms_, basis_);
268 const auto nao = BasisSet::nbf(basis_);
270 assert(nao == coefs_.rows());
271 const auto shell2ao = BasisSet::compute_shell2bf(basis_);
273 for (
size_t iatom = 0; iatom < atoms_.size(); ++iatom) {
274 for (
auto ish : atom2shell_[iatom]) {
275 auto ao = shell2ao[ish];
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;
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;
287 END_FOR_SOLIDHARM_MOLDEN
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;
333 template <
typename Coeffs,
typename Occs,
typename Energies = Eigen::VectorXd>
335 const Eigen::VectorXd& cell_axes,
336 const std::vector<Shell>& basis,
337 const Coeffs& coefficients,
338 const Occs& occupancies,
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)
353 os <<
"[SpaceGroup] (Number)" << std::endl;
354 os << space_group_ << std::endl;
359 os <<
"[Operators]" << std::endl;
360 os <<
"x, y, z" << std::endl;
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;
372 void write(
const std::string& filename)
const {
373 std::ofstream os(filename);
384 Eigen::VectorXd cell_axes_;
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