Skip to content

Commit

Permalink
Modernize BasisSet::print_detail_cfour (#2937)
Browse files Browse the repository at this point in the history 
* clang-format

* update docstrings and URL

* make string const

* int64_t-ifying pass

* remove sprintf usage, part 1

* add none_of_equal utility function, rework complicated std::find expression

* avoid using == for floating point values

* finished

* hopefully fix build error, add docstring, restrict none_of_equal from processing non-double FP types

* capture 'this'

* add missing include

* the two new fns are probably better off not being member fns

* also update rotten link in the native python implementation

* add one more lambda to make an object const

* the decltype magic was probably excessive
  • Loading branch information
@TiborGY
TiborGY committed Nov 14, 2023
1 parent 289d095 commit 0d54746
Show file tree
Hide file tree
Showing 3 changed files with 91 additions and 73 deletions.
2 changes: 1 addition & 1 deletion psi4/driver/qcdb/libmintsbasisset.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 = ''
Expand Down
121 changes: 54 additions & 67 deletions psi4/src/psi4/libmints/basisset.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@
#include <regex>
#include <stdexcept>
#include <cstdio>
#include <iomanip>
#include <cstdlib>
#include <cmath>
#include <map>
Expand Down Expand Up @@ -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<std::vector<int>> 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<std::vector<int>> 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<std::vector<double>> exp_per_am(max_am_center + 1);
std::vector<std::vector<double>> 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++) {
Expand All @@ -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();
Expand Down Expand Up @@ -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!");
Expand Down
41 changes: 36 additions & 5 deletions psi4/src/psi4/libmints/basisset.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@
#include <string>
#include <vector>
#include <map>
#include <type_traits>
#include <algorithm>
PRAGMA_WARNING_PUSH
PRAGMA_WARNING_IGNORE_DEPRECATED_DECLARATIONS
#include <memory>
Expand Down Expand Up @@ -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_.
Expand Down Expand Up @@ -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 <typename T>
bool none_of_equal(const std::vector<T> &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<T>::value && !std::is_same<T, double>::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<T>::value) {
return fpeq(X, value);
} else {
return X == value;
}
});
}

} // namespace psi

#endif

0 comments on commit 0d54746

Please sign in to comment.