evaleev / libint

Libint: high-performance library for computing Gaussian integrals in quantum mechanics
Other
229 stars 95 forks source link

Buffer size of Engine Operator::nuclear derivatives appears to be too large #199

Closed adabbott closed 3 years ago

adabbott commented 3 years ago

It appears as though the buffer size when calling first and second derivatives of Operator::nuclear integrals is larger than it needs to be. For example, for H2 nuclear-electron attraction integrals, the buffer is size 18 when there are only 12 derivatives, 6 w.r.t. shell coordinates and 6 w.r.t. nuclear coordinates in the operator. The below code snippet shows this behavior:

Click to expand code snippet

```c++ #include #include // Computes nuclear derivatives of potential energy integrals int main(int argc, char* argv[]) { using namespace libint2; using namespace std; libint2::initialize(); string xyzfilename = "/home/adabbott/h2.xyz"; ifstream input_file(xyzfilename); vector atoms = read_dotxyz(input_file); BasisSet obs("sto-3g", atoms); int nbf = obs.nbf(); int deriv_order = 1; // Nuclear attraction integral derivative engine Engine v_engine(Operator::nuclear,obs.max_nprim(),obs.max_l(),deriv_order); v_engine.set_params(make_point_charges(atoms)); const auto& buf_vec = v_engine.results(); // Print buffer vector size and theoretical buffer vector size according to libint wiki // For k centers and l atoms // (3 * (k + l) + 0) * (3 * (k + l) + 1) * ... * (3 * (k+l) + n - 1) / n! int theo_buf_vec_size = 1; int factorial = 1; for (auto i=0; i

Generally, for first derivatives, the leading index of the buffer appears to be 6 + NCART + NCART, for NCART cartesian coordinates in the operator. The first NCART is empty, just zeros, and the last NCART is the true operator derivatives. This doesn't appear to be intended, given this comment in engine.impl.h and the wording in the wiki.

For second derivatives, the leading buffer index is the flattened upper triangle of an (6 + NCART + NCART) by (6 + NCART + NCART) array, leading to 171 elements in the case of H2 when only 78 should be needed. For larger molecules, this results in buffers that are unnecessarily huge.

Is there a reason for this I am not catching?