diff --git a/psi4/driver/qcdb/libmintsbasisset.py b/psi4/driver/qcdb/libmintsbasisset.py index 9091bd3fd5b..0b564e7d9fd 100644 --- a/psi4/driver/qcdb/libmintsbasisset.py +++ b/psi4/driver/qcdb/libmintsbasisset.py @@ -1360,7 +1360,7 @@ def print_detail_gamess(self, out=None, numbersonly=False): def print_detail_cfour(self, out=None): """Returns a string in CFOUR-style of the basis (per-atom) - * Format from http://slater.chemie.uni-mainz.de/cfour/index.php?n=Main.OldFormatOfAnEntryInTheGENBASFile + * Format from https://web.archive.org/web/20221130013041/http://slater.chemie.uni-mainz.de/cfour/index.php?n=Main.OldFormatOfAnEntryInTheGENBASFile """ text = '' diff --git a/psi4/src/psi4/libmints/basisset.cc b/psi4/src/psi4/libmints/basisset.cc index 9d2f4204125..ec2c9efaf14 100644 --- a/psi4/src/psi4/libmints/basisset.cc +++ b/psi4/src/psi4/libmints/basisset.cc @@ -56,6 +56,7 @@ #include #include #include +#include #include #include #include @@ -476,66 +477,55 @@ void BasisSet::print_detail(std::string out) const { } } +/// @brief Returns a string in CFOUR-style of the basis (per-atom). Format from +/// https://web.archive.org/web/20221130013041/http://slater.chemie.uni-mainz.de/cfour/index.php?n=Main.OldFormatOfAnEntryInTheGENBASFile +/// @return CFOUR-style of the basis (per-atom) std::string BasisSet::print_detail_cfour() const { - char buffer[120]; std::stringstream ss; - std::string nameUpperCase = name_; - to_upper(nameUpperCase); + ss << std::fixed << std::showpoint; + const std::string nameUpperCase = to_upper_copy(name_); for (int uA = 0; uA < molecule_->nunique(); uA++) { - const int A = molecule_->unique(uA); - - sprintf(buffer, "%s:P4_%d\n", molecule_->symbol(A).c_str(), A + 1); - ss << buffer; - sprintf(buffer, "Psi4 basis %s for element %s atom %d\n\n", nameUpperCase.c_str(), molecule_->symbol(A).c_str(), - A + 1); - ss << buffer; - - int first_shell = center_to_shell_[A]; - int n_shell = center_to_nshell_[A]; - - int max_am_center = 0; - for (int Q = 0; Q < n_shell; Q++) - max_am_center = - (shells_[Q + first_shell].am() > max_am_center) ? shells_[Q + first_shell].am() : max_am_center; - - std::vector> shell_per_am(max_am_center + 1); - for (int Q = 0; Q < n_shell; Q++) shell_per_am[shells_[Q + first_shell].am()].push_back(Q); + const auto A = molecule_->unique(uA); + ss << molecule_->symbol(A) << ":P4_" << A + 1 << '\n'; + ss << "Psi4 basis " << nameUpperCase << " for element " << molecule_->symbol(A) << " atom " << A + 1 << "\n\n"; + + const auto first_shell = center_to_shell_[A]; + const auto n_shell = center_to_nshell_[A]; + // Use immediately evaluated lambdas for complicated initialiation of const objects + const auto max_am_center = [&] { + int max = 0; + for (int Q = 0; Q < n_shell; Q++) { + if (shells_[Q + first_shell].am() > max) max = shells_[Q + first_shell].am(); + } + return max; + }(); + const auto shell_per_am = [&] { + std::vector> tmpvec(max_am_center + 1); + for (int Q = 0; Q < n_shell; Q++) tmpvec[shells_[Q + first_shell].am()].push_back(Q); + return tmpvec; + }(); // Write number of shells in the basis set - sprintf(buffer, "%3d\n", max_am_center + 1); - ss << buffer; - + ss << to_str_width(max_am_center + 1, 3) << '\n'; // Write angular momentum for each shell - for (int am = 0; am <= max_am_center; am++) { - sprintf(buffer, "%5d", am); - ss << buffer; - } - sprintf(buffer, "\n"); - ss << buffer; - + for (int am = 0; am <= max_am_center; am++) ss << to_str_width(am, 5); + ss << '\n'; // Write number of contracted basis functions for each shell - for (int am = 0; am <= max_am_center; am++) { - sprintf(buffer, "%5lu", shell_per_am[am].size()); - ss << buffer; - } - sprintf(buffer, "\n"); - ss << buffer; + for (int am = 0; am <= max_am_center; am++) ss << to_str_width(shell_per_am[am].size(), 5); + ss << '\n'; std::vector> exp_per_am(max_am_center + 1); std::vector> coef_per_am(max_am_center + 1); for (int am = 0; am <= max_am_center; am++) { - // TODO: std::find safe on floats? seems to work // Collect unique exponents among all functions for (size_t Q = 0; Q < shell_per_am[am].size(); Q++) { for (int K = 0; K < shells_[shell_per_am[am][Q] + first_shell].nprimitive(); K++) { - if (!(std::find(exp_per_am[am].begin(), exp_per_am[am].end(), - shells_[shell_per_am[am][Q] + first_shell].exp(K)) != exp_per_am[am].end())) { + if (none_of_equal(exp_per_am[am], shells_[shell_per_am[am][Q] + first_shell].exp(K))) { exp_per_am[am].push_back(shells_[shell_per_am[am][Q] + first_shell].exp(K)); } } } - // Collect coefficients for each exp among all functions, zero otherwise for (size_t Q = 0; Q < shell_per_am[am].size(); Q++) { for (size_t ep = 0, K = 0; ep < exp_per_am[am].size(); ep++) { @@ -552,44 +542,33 @@ std::string BasisSet::print_detail_cfour() const { } // Write number of exponents for each shell - for (int am = 0; am <= max_am_center; am++) { - sprintf(buffer, "%5lu", exp_per_am[am].size()); - ss << buffer; - } - sprintf(buffer, "\n\n"); - ss << buffer; + for (int am = 0; am <= max_am_center; am++) ss << to_str_width(exp_per_am[am].size(), 5); + ss << "\n\n"; for (int am = 0; am <= max_am_center; am++) { // Write exponents for each shell for (size_t ep = 0; ep < exp_per_am[am].size(); ep++) { - if (exp_per_am[am][ep] >= 10000000.0) - sprintf(buffer, "%13.4f ", exp_per_am[am][ep]); - else if (exp_per_am[am][ep] >= 1000000.0) - sprintf(buffer, "%13.5f ", exp_per_am[am][ep]); - else if (exp_per_am[am][ep] >= 100000.0) - sprintf(buffer, "%13.6f ", exp_per_am[am][ep]); - else - sprintf(buffer, "%14.7f", exp_per_am[am][ep]); - ss << buffer; - if (((ep + 1) % 5 == 0) || ((ep + 1) == exp_per_am[am].size())) { - sprintf(buffer, "\n"); - ss << buffer; + if (exp_per_am[am][ep] >= 10000000.0) { + ss << std::setprecision(4) << std::setw(13) << exp_per_am[am][ep] << ' '; + } else if (exp_per_am[am][ep] >= 1000000.0) { + ss << std::setprecision(5) << std::setw(13) << exp_per_am[am][ep] << ' '; + } else if (exp_per_am[am][ep] >= 100000.0) { + ss << std::setprecision(6) << std::setw(13) << exp_per_am[am][ep] << ' '; + } else { + ss << std::setprecision(7) << std::setw(14) << exp_per_am[am][ep] << ' '; } + if (((ep + 1) % 5 == 0) || ((ep + 1) == exp_per_am[am].size())) ss << '\n'; } - sprintf(buffer, "\n"); - ss << buffer; + ss << '\n'; // Write contraction coefficients for each shell for (size_t ep = 0; ep < exp_per_am[am].size(); ep++) { for (size_t bf = 0; bf < shell_per_am[am].size(); bf++) { - sprintf(buffer, "%10.7f ", coef_per_am[am][bf * exp_per_am[am].size() + ep]); - ss << buffer; + ss << std::setprecision(7) << std::setw(10) << coef_per_am[am][bf * exp_per_am[am].size() + ep] << ' '; } - sprintf(buffer, "\n"); - ss << buffer; + ss << '\n'; } - sprintf(buffer, "\n"); - ss << buffer; + ss << '\n'; } } return ss.str(); @@ -1273,6 +1252,14 @@ void BasisSet::compute_phi(double *phi_ao, double x, double y, double z) { } // nshell } +bool psi::fpeq(const double a, const double b, const double THR/* = 1.0E-14*/) { + if (std::abs(a - b) < THR) { + return true; + } else { + return false; + } +} + void BasisSet::convert_sap_contraction() { if(max_am_ != 0) { throw PSIEXCEPTION("SAP potentials should be composed of a single S function per atom, and not contain higher angular momentum!"); diff --git a/psi4/src/psi4/libmints/basisset.h b/psi4/src/psi4/libmints/basisset.h index ef544a319a1..bf8a984d1e8 100644 --- a/psi4/src/psi4/libmints/basisset.h +++ b/psi4/src/psi4/libmints/basisset.h @@ -42,6 +42,8 @@ #include #include #include +#include +#include PRAGMA_WARNING_PUSH PRAGMA_WARNING_IGNORE_DEPRECATED_DECLARATIONS #include @@ -346,9 +348,9 @@ class PSI_API BasisSet { void print_detail(std::string out) const; void print_detail() const { print_detail("outfile"); } - /** Returns a string in CFOUR-style of the basis (per-atom) - * Format from http://slater.chemie.uni-mainz.de/cfour/index.php?n=Main.OldFormatOfAnEntryInTheGENBASFile - */ + /// @brief Returns a string in CFOUR-style of the basis (per-atom). Format from + /// https://web.archive.org/web/20221130013041/http://slater.chemie.uni-mainz.de/cfour/index.php?n=Main.OldFormatOfAnEntryInTheGENBASFile + /// @return CFOUR-style of the basis (per-atom) std::string print_detail_cfour() const; /** Refresh internal basis set data. Useful if someone has pushed to shells_. @@ -405,16 +407,45 @@ class PSI_API BasisSet { void move_atom(int atom, const Vector3 &trans); // Returns the values of the basis functions at a point void compute_phi(double *phi_ao, double x, double y, double z); - // Converts the contraction to match the SAP approach. void convert_sap_contraction(); - private: + private: /// Helper functions for frozen core to reduce LOC int atom_to_period(int Z); int period_to_full_shell(int p); }; +/// @brief Determine if two floating-point numbers are practically equal. Simply comparing two FP numbers is usually +/// a bad idea due to numerical noise. They *can* be equal, but only in fairly specific circumstances. +/// @param a : first number +/// @param b : second number +/// @param THR : (optional) equality threshold +/// @return True if their absolute difference is smaller than THR, false otherwise. +bool fpeq(const double a, const double b, const double THR = 1.0E-14); + +/// @brief Determine if a particular value is abscent from an std::vector. If the type is floating-point, make sure +/// it is a double and use a fuzzy equality to account for numerical noise. +/// @tparam T : Type of the value. If floating-point, it must be a double. +/// @param container : Containter to search in +/// @param value : Value to look for +/// @return : true if none of the elements of container are equal to value +template +bool none_of_equal(const std::vector &container, const T value) { + // The equality threshold in fpeq is set for the expected precision of double. Other FP types could be + // implemented in the future if needed, but this function should refuse to process them until then. + // Without this check implicit conversion rules could silently promote eg. a float to double when calling fpeq. + static_assert(false == (std::is_floating_point::value && !std::is_same::value), + "Support for FP types other than double is not implemented in none_of_equal"); + return std::none_of(container.cbegin(), container.cend(), [value](const T X) { + if constexpr (std::is_floating_point::value) { + return fpeq(X, value); + } else { + return X == value; + } + }); +} + } // namespace psi #endif