• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

adc-connect / adcc / 11687183074

05 Nov 2024 03:20PM UTC coverage: 73.291% (-0.5%) from 73.743%
11687183074

Pull #169

github

web-flow
Merge e14f10df3 into e2e1bce49
Pull Request #169: Additional operator integrals and selection of gauge origin

1166 of 1868 branches covered (62.42%)

Branch coverage included in aggregate %.

132 of 233 new or added lines in 12 files covered. (56.65%)

5 existing lines in 3 files now uncovered.

7313 of 9701 relevant lines covered (75.38%)

271899.11 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.69
/libadcc/ReferenceState.cc
1
//
2
// Copyright (C) 2019 by the adcc authors
3
//
4
// This file is part of adcc.
5
//
6
// adcc is free software: you can redistribute it and/or modify
7
// it under the terms of the GNU General Public License as published
8
// by the Free Software Foundation, either version 3 of the License, or
9
// (at your option) any later version.
10
//
11
// adcc 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 General Public License for more details.
15
//
16
// You should have received a copy of the GNU General Public License
17
// along with adcc. If not, see <http://www.gnu.org/licenses/>.
18
//
19

20
#include "ReferenceState.hh"
21
#include "MoIndexTranslation.hh"
22
#include "TensorImpl.hh"
23
#include "exceptions.hh"
24
#include "import_eri.hh"
25
#include "make_symmetry.hh"
26

27
namespace libadcc {
28
ReferenceState::ReferenceState(std::shared_ptr<const HartreeFockSolution_i> hfsoln_ptr,
668✔
29
                               std::shared_ptr<const MoSpaces> mo_ptr,
30
                               bool symmetry_check_on_import)
668✔
31
      : m_n_alpha{0},
668✔
32
        m_n_beta{0},
668✔
33
        m_orben{},
668✔
34
        m_orbcoeff{},
668✔
35
        m_orbcoeff_alpha{},
668✔
36
        m_orbcoeff_beta{},
668✔
37
        m_fock{},
668✔
38
        m_eri{},
668✔
39
        m_hfsoln_ptr{hfsoln_ptr},
668✔
40
        m_mo_ptr{mo_ptr},
668✔
41
        m_symmetry_check_on_import{symmetry_check_on_import},
668✔
42
        m_timer{} {
668✔
43
  for (auto& ss : mo_ptr->subspaces_occupied) {
1,702✔
44
    m_n_alpha += mo_ptr->n_orbs_alpha(ss);
1,034✔
45
    m_n_beta += mo_ptr->n_orbs_beta(ss);
1,034✔
46
  }
47

48
  // Import orbital energies and coefficients
49
  import_orbital_energies(*hfsoln_ptr, mo_ptr, symmetry_check_on_import);
668✔
50
  import_orbital_coefficients(*hfsoln_ptr, mo_ptr, symmetry_check_on_import);
668✔
51
}
668✔
52

53
//
54
// Importers
55
//
56
void ReferenceState::import_orbital_energies(
668✔
57
      const HartreeFockSolution_i& hf, const std::shared_ptr<const MoSpaces>& mo_ptr,
58
      bool symmetry_check_on_import) {
59
  m_orben.clear();
668✔
60
  RecordTime rec(m_timer, "import/orbital_energies");
1,336✔
61

62
  std::vector<scalar_type> tmp_orben(hf.n_orbs());
668✔
63
  hf.orben_f(tmp_orben.data(), tmp_orben.size());
668✔
64
  for (auto& ss : mo_ptr->subspaces) {
2,502✔
65
    std::shared_ptr<Symmetry> sym_ptr = make_symmetry_orbital_energies(mo_ptr, ss);
1,834✔
66
    std::shared_ptr<Tensor> orben     = make_tensor_zero(sym_ptr);
1,834✔
67
    MoIndexTranslation idxtrans(mo_ptr, ss);
1,834✔
68

69
    auto orben_generator = [&tmp_orben, &idxtrans](
2,362✔
70
                                 const std::vector<std::pair<size_t, size_t>>& block,
71
                                 scalar_type* buffer) {
7,094✔
72
      const std::pair<size_t, size_t>& block_ss = block[0];
2,362✔
73
      const size_t length_ss_full               = block_ss.second - block_ss.first;
2,362✔
74

75
      std::vector<RangeMapping> mappings = idxtrans.map_range_to_hf_provider(block);
4,724✔
76
      for (const RangeMapping& map : mappings) {
4,728✔
77
        // orben is 1D -> idim == 0 always
78
        const size_t idim = 0;
2,366✔
79

80
        // Starting index of the mapping in the adcc index convention
81
        const size_t start = map.from().axis(idim).start() - block_ss.first;
2,366✔
82

83
        // Start and end in the HfProvider indexing convention
84
        const size_t start_hf = map.to().axis(idim).start();
2,366✔
85
        const size_t end_hf   = map.to().axis(idim).end();
2,366✔
86

87
        if (start + (end_hf - start_hf) > length_ss_full) {
2,366✔
88
          throw runtime_error("Internal error: Out-of-bounds write in orben_generator");
×
89
        }
90
        // Do the "import" by copying the appropriate part of the buf_orben object
91
        std::copy(tmp_orben.begin() + static_cast<ptrdiff_t>(start_hf),
2,366✔
92
                  tmp_orben.begin() + static_cast<ptrdiff_t>(end_hf), buffer + start);
4,732✔
93
      }
94
    };
2,362✔
95
    orben->import_from(orben_generator, hf.conv_tol(), symmetry_check_on_import);
1,834✔
96
    orben->set_immutable();
1,834✔
97
    m_orben[ss] = orben;
1,834✔
98
  }
1,834✔
99
}
668✔
100

101
void ReferenceState::import_orbital_coefficients(
668✔
102
      const HartreeFockSolution_i& hf, const std::shared_ptr<const MoSpaces>& mo_ptr,
103
      bool symmetry_check_on_import) {
104
  RecordTime rec(m_timer, "import/orbital_coefficients");
1,336✔
105
  m_orbcoeff.clear();
668✔
106
  m_orbcoeff_alpha.clear();
668✔
107
  m_orbcoeff_beta.clear();
668✔
108

109
  std::vector<scalar_type> tmp_orbcoeff(hf.n_orbs() * hf.n_bas());
668✔
110
  tmp_orbcoeff.resize(hf.n_orbs() * hf.n_bas());
668✔
111
  hf.orbcoeff_fb(tmp_orbcoeff.data(), tmp_orbcoeff.size());
668✔
112
  for (auto& ss : mo_ptr->subspaces) {
2,502✔
113
    // The index translation for the first subspace (which is not "b")
114
    MoIndexTranslation idxtrans(mo_ptr, ss);
1,834✔
115

116
    auto orbcoeff_generator_ab = [&tmp_orbcoeff, &idxtrans, &hf](
5,772✔
117
                                       const std::vector<std::pair<size_t, size_t>>&
118
                                             block,
119
                                       scalar_type* buffer, bool alpha, bool beta) {
52,012✔
120
      const std::pair<size_t, size_t>& block_ss = block[0];
5,772✔
121
      const std::pair<size_t, size_t>& block_b  = block[1];
5,772✔
122

123
      // Length of the axis:
124
      const size_t length_b       = block_b.second - block_b.first;
5,772✔
125
      const size_t length_ss_full = block_ss.second - block_ss.first;
5,772✔
126

127
      std::vector<RangeMapping> mappings_ss =
128
            idxtrans.map_range_to_hf_provider({{block_ss}});
11,544✔
129
      for (const RangeMapping& map_ss : mappings_ss) {
11,552✔
130
        // Get start and end of index ranges in adcc and HfProvider convention
131
        const size_t idim      = 0;  // The ss subspace of orbcoeff is 1D => idim == 0
5,780✔
132
        const size_t start_ss  = map_ss.from().axis(idim).start();
5,780✔
133
        const size_t end_ss    = map_ss.from().axis(idim).end();
5,780✔
134
        const size_t length_ss = end_ss - start_ss;
5,780✔
135

136
        const size_t start_hf = map_ss.to().axis(idim).start();
5,780✔
137
        const size_t end_hf   = map_ss.to().axis(idim).end();
5,780✔
138

139
        // Since the orbecoeff are 2D-arrays and moreover the tmp_orbcoeff vector is a 2D
140
        // array in row-major C-like index ordering, the first element of this buffer to
141
        // copy is
142
        auto ittmp = tmp_orbcoeff.begin() +
11,560✔
143
                     static_cast<ptrdiff_t>(start_hf * hf.n_bas() + block_b.first);
5,780✔
144

145
        // The first element into which data will be copied, on the other hand,
146
        // is shifted, like the following, since in the output buffer there will
147
        // length_b columns and the first (start_ss - block_ss.first) rows are skipped
148
        // in this mapping loop.
149
        scalar_type* ptrbuf =
5,780✔
150
              buffer + static_cast<ptrdiff_t>((start_ss - block_ss.first) * length_b);
5,780✔
151

152
        if (alpha && beta) {
5,780✔
153
          throw runtime_error("Internal error: Requested both alpha and beta block.");
×
154
        }
155

156
        // Fill with zeros if we are in the alpha or beta block and this
157
        // should be skipped
158
        bool fill_zeros = false;
5,780✔
159
        if (idxtrans.spin_of({start_ss}) == "b") {
5,780✔
160
          if (idxtrans.spin_of({end_ss - 1}) != "b") {
2,894✔
161
            throw runtime_error(
×
162
                  "Internal error: Assertion about blocks being split whenever spin "
163
                  "changes is violated.");
×
164
          }
165
          if (!beta) {
2,894✔
166
            // beta block should be zero
167
            fill_zeros = true;
582✔
168
          }
169
        }
170
        if (idxtrans.spin_of({end_ss - 1}) == "a") {
5,780✔
171
          if (idxtrans.spin_of({start_ss}) != "a") {
2,886✔
172
            throw runtime_error(
×
173
                  "Internal error: Assertion about blocks being split whenever spin "
174
                  "changes is violated.");
×
175
          }
176
          if (!alpha) {
2,886✔
177
            // beta block should be zero
178
            fill_zeros = true;
578✔
179
          }
180
        }
181

182
        const ptrdiff_t slength_b = static_cast<ptrdiff_t>(length_b);
5,780✔
183
        const ptrdiff_t sn_bas    = static_cast<ptrdiff_t>(hf.n_bas());
5,780✔
184

185
        // Copy length_b elements at a time (one row) and then advance the pointer into
186
        // the buffer (ptrbuf) and the iterator into the tmp_orbcoeff vector.
187
        // For the ittmp iterator, n_bas() elements need to be skipped, for ptrbuf only
188
        // length_b, because it only is a matrix block with length_b columns.
189
        for (size_t irow = 0; irow < length_ss;
34,244✔
190
             ++irow, ptrbuf += length_b, ittmp += sn_bas) {
28,464✔
191
          if (ptrbuf + slength_b > buffer + length_ss_full * length_b) {
28,464✔
192
            throw runtime_error(
×
193
                  "Internal error: Out-of-bounds write in orbcoeff_generator");
×
194
          }
195
          if (fill_zeros) {
28,464✔
196
            std::fill(ptrbuf, ptrbuf + length_b, 0);
5,892✔
197
          } else {
198
            std::copy(ittmp, ittmp + slength_b, ptrbuf);
22,572✔
199
          }
200
        }
201

202
        if (ittmp != tmp_orbcoeff.begin() +
5,780✔
203
                           static_cast<ptrdiff_t>(end_hf * hf.n_bas() + block_b.first)) {
5,780✔
204
          throw runtime_error("Internal error: ittmp index issue in orbcoeff_generator");
×
205
        }
206
        if (ptrbuf != buffer + (end_ss - block_ss.first) * length_b) {
5,780✔
207
          throw runtime_error("Internal error: ptrbuf index issue in orbcoeff_generator");
×
208
        }
209
      }
210
    };
5,772✔
211

212
    //
213
    // Import alpha coefficients, i.e. the vector (a)
214
    //                                            (0)
215
    //
216
    std::shared_ptr<Symmetry> sym_a_ptr =
217
          make_symmetry_orbital_coefficients(mo_ptr, ss + "b", hf.n_bas(), "a");
3,668✔
218
    std::shared_ptr<Tensor> orbcoeff_a = make_tensor_zero(sym_a_ptr);
1,834✔
219
    orbcoeff_a->import_from(
3,668✔
220
          [&orbcoeff_generator_ab](const std::vector<std::pair<size_t, size_t>>& block,
2,886✔
221
                                   scalar_type* buffer) {
2,886✔
222
            orbcoeff_generator_ab(block, buffer, /*alpha*/ true, /*beta*/ false);
2,886✔
223
          },
2,886✔
224
          hf.conv_tol(), symmetry_check_on_import);
1,834✔
225
    orbcoeff_a->set_immutable();
1,834✔
226
    m_orbcoeff_alpha[ss + "b"] = orbcoeff_a;
1,834✔
227

228
    //
229
    // Import beta coefficients, i.e. the vector (0)
230
    //                                           (b)
231
    //
232
    std::shared_ptr<Symmetry> sym_b_ptr =
233
          make_symmetry_orbital_coefficients(mo_ptr, ss + "b", hf.n_bas(), "b");
3,668✔
234
    std::shared_ptr<Tensor> orbcoeff_b = make_tensor_zero(sym_b_ptr);
1,834✔
235
    orbcoeff_b->import_from(
3,668✔
236
          [&orbcoeff_generator_ab](const std::vector<std::pair<size_t, size_t>>& block,
2,886✔
237
                                   scalar_type* buffer) {
2,886✔
238
            orbcoeff_generator_ab(block, buffer, /*alpha*/ false, /*beta*/ true);
2,886✔
239
          },
2,886✔
240
          hf.conv_tol(), symmetry_check_on_import);
1,834✔
241
    orbcoeff_b->set_immutable();
1,834✔
242
    m_orbcoeff_beta[ss + "b"] = orbcoeff_b;
1,834✔
243

244
    //
245
    // Build full coefficients (a 0)
246
    //                         (0 b)
247
    //
248
    // Export orbcoeff_a and orbcoeff_b for the full coefficient import
249
    // (The rationale is that these already have the right MO ordering along the MO axis,
250
    //  which makes the code for the full coefficient import a lot more readable)
251
    std::vector<scalar_type> exp_orbcoeff_a;
1,834✔
252
    orbcoeff_a->export_to(exp_orbcoeff_a);
1,834✔
253
    std::vector<scalar_type> exp_orbcoeff_b;
1,834✔
254
    orbcoeff_b->export_to(exp_orbcoeff_b);
1,834✔
255
    if (exp_orbcoeff_a.size() != exp_orbcoeff_b.size()) {
1,834✔
256
      throw runtime_error(
×
257
            "Internal error: exp_orbcoeff_a.size() and exp_orbcoeff_b.size() do not "
258
            "agree.");
×
259
    }
260

261
    auto orbcoeff_generator = [&exp_orbcoeff_a, &exp_orbcoeff_b, &hf](
4,042✔
262
                                    const std::vector<std::pair<size_t, size_t>>& block,
263
                                    scalar_type* buffer) {
72,584✔
264
      const std::pair<size_t, size_t>& block_ss = block[0];
4,042✔
265
      const std::pair<size_t, size_t>& block_b  = block[1];
4,042✔
266

267
      // Length of the axis:
268
      const size_t length_b       = block_b.second - block_b.first;
4,042✔
269
      const size_t length_ss_full = block_ss.second - block_ss.first;
4,042✔
270

271
      // Spin data for axis b
272
      const std::vector<std::pair<size_t, size_t>> b_spin_ranges{
273
            {0, hf.n_bas()}, {hf.n_bas(), 2 * hf.n_bas()}};
4,042✔
274
      const std::vector<scalar_type*> b_spin_data{exp_orbcoeff_a.data(),
4,042✔
275
                                                  exp_orbcoeff_b.data()};
4,042✔
276

277
      // Loop over spins (alpha == 0 and beta == 1)
278
      for (size_t ispin = 0; ispin < 2; ++ispin) {
12,126✔
279
        size_t from_b_first  = block_b.first;
8,084✔
280
        size_t from_b_length = length_b;
8,084✔
281

282
        if (block_b.first >= b_spin_ranges[ispin].second) {
8,084✔
283
          // Beyond the range covered by this spin
284
          continue;
1,156✔
285
        }
286
        if (block_b.second - 1 < b_spin_ranges[ispin].first) {
6,928✔
287
          // Before the range covered by this spin
288
          continue;
2,886✔
289
        }
290

291
        // Get the first index in the b axis to be used and the length
292
        // for the spin chunk covered in this loop iteration
293
        from_b_first  = std::max(block_b.first, b_spin_ranges[ispin].first);
4,042✔
294
        from_b_length = std::min(hf.n_bas(), block_b.second - from_b_first);
4,042✔
295

296
        // from_b_first has still spin information (i.e. indexed in the range
297
        // [0,2*n_bas], but exp_orbcoeff_a and exp_orbcoeff_b containeres
298
        // only covers the spatial range [0, n_bas], so we reduce it to
299
        // the spatial index only
300
        const size_t from_b_first_spatial = from_b_first % hf.n_bas();
4,042✔
301

302
        // Prepare input and output pointers and copy the data.
303
        scalar_type* from   = b_spin_data[ispin];
4,042✔
304
        scalar_type* ptrbuf = buffer + (from_b_first - block_b.first) * length_ss_full;
4,042✔
305
        size_t offset       = hf.n_bas() * block_ss.first + from_b_first_spatial;
4,042✔
306
        for (size_t iss = 0; iss < length_ss_full;
24,166✔
307
             ++iss, ptrbuf += length_b, offset += hf.n_bas()) {
20,124✔
308
          if (offset + from_b_length > exp_orbcoeff_a.size()) {
20,124✔
309
            throw runtime_error("Read buffer overflow.");
×
310
          }
311
          if (ptrbuf + from_b_length > buffer + length_b * length_ss_full) {
20,124✔
312
            throw runtime_error("Write buffer overflow.");
×
313
          }
314
          if (ptrbuf < buffer) {
20,124✔
315
            throw runtime_error("Write buffer underflow.");
×
316
          }
317
          std::copy(from + offset, from + offset + from_b_length, ptrbuf);
20,124✔
318
        }
319
      }
320
    };
4,042✔
321
    std::shared_ptr<Symmetry> sym_ab_ptr =
322
          make_symmetry_orbital_coefficients(mo_ptr, ss + "b", hf.n_bas(), "ab");
3,668✔
323
    std::shared_ptr<Tensor> orbcoeff_ab = make_tensor_zero(sym_ab_ptr);
1,834✔
324
    orbcoeff_ab->import_from(orbcoeff_generator, hf.conv_tol(), symmetry_check_on_import);
1,834✔
325
    orbcoeff_ab->set_immutable();
1,834✔
326
    m_orbcoeff[ss + "b"] = orbcoeff_ab;
1,834✔
327
  }
1,834✔
328
}
668✔
329

330
//
331
// Info about the reference state
332
//
333
size_t ReferenceState::spin_multiplicity() const {
98✔
334
  if (!restricted()) {
98✔
335
    return 0;
48✔
336
  } else {
337
    return m_hfsoln_ptr->spin_multiplicity();
50✔
338
  }
339
}
340

341
std::string ReferenceState::irreducible_representation() const {
98✔
342
  if (m_mo_ptr->point_group == "C1") return "A";
98✔
343

344
  // Notice: This is not true, since the occupation might be uneven
345
  // between spatial parts
346
  // if (n_alpha() == n_beta()) return m_mo_ptr->irrep_totsym();
347

348
  // For non-closed-shell we would need to see whatever orbitals are occupied
349
  // and determine from this the resulting overall irrep
350
  throw not_implemented_error("Only C1 is implemented.");
×
351
}
352

353
std::vector<scalar_type> ReferenceState::nuclear_multipole(
462✔
354
      size_t order, std::vector<scalar_type> gauge_origin) const {
355
  std::vector<scalar_type> ret((order + 2) * (order + 1) / 2);
462✔
356
  m_hfsoln_ptr->nuclear_multipole(order, gauge_origin, ret.data(), ret.size());
462✔
357
  return ret;
462✔
NEW
358
}
×
359

NEW
360
std::vector<scalar_type> ReferenceState::determine_gauge_origin(
×
361
      std::string gauge_origin) const {
NEW
362
  std::vector<scalar_type> ret = m_hfsoln_ptr->determine_gauge_origin(gauge_origin);
×
363
  return ret;
×
364
}
365

366
//
367
// Tensor data
368
//
369
std::shared_ptr<Tensor> ReferenceState::orbital_energies(const std::string& space) const {
326✔
370
  const auto itfound = m_orben.find(space);
326✔
371
  if (itfound == m_orben.end()) {
326✔
372
    throw invalid_argument(
×
373
          "Invalid space string " + space + ": An object orbital_energies(" + space +
×
374
          ") is not known. Check that the string has exactly 2 characters "
375
          "and describes a valid orbital subspace.");
×
376
  }
377
  return itfound->second;
652✔
378
}
379

380
std::shared_ptr<Tensor> ReferenceState::orbital_coefficients(
130,654✔
381
      const std::string& space) const {
382
  const auto itfound = m_orbcoeff.find(space);
130,654✔
383
  if (itfound == m_orbcoeff.end()) {
130,654✔
384
    throw invalid_argument(
×
385
          "Invalid space string " + space + ": An object orbital_coefficients(" + space +
×
386
          ") is not known. Check that the string has exactly 3 characters "
387
          "and describes two valid orbital subspace.");
×
388
  }
389
  return itfound->second;
261,308✔
390
}
391

392
std::shared_ptr<Tensor> ReferenceState::orbital_coefficients_alpha(
2,576✔
393
      const std::string& space) const {
394
  const auto itfound = m_orbcoeff_alpha.find(space);
2,576✔
395
  if (itfound == m_orbcoeff_alpha.end()) {
2,576✔
396
    throw invalid_argument(
×
397
          "Invalid space string " + space + ": An object orbital_coefficients_alpha(" +
×
398
          space +
×
399
          ") is not known. Check that the string has exactly 3 characters "
400
          "and describes two valid orbital subspace.");
×
401
  }
402
  return itfound->second;
5,152✔
403
}
404

405
std::shared_ptr<Tensor> ReferenceState::orbital_coefficients_beta(
2,576✔
406
      const std::string& space) const {
407
  const auto itfound = m_orbcoeff_beta.find(space);
2,576✔
408
  if (itfound == m_orbcoeff_beta.end()) {
2,576✔
409
    throw invalid_argument(
×
410
          "Invalid space string " + space + ": An object orbital_coefficients_beta(" +
×
411
          space +
×
412
          ") is not known. Check that the string has exactly 3 characters "
413
          "and describes two valid orbital subspace.");
×
414
  }
415
  return itfound->second;
5,152✔
416
}
417

418
std::shared_ptr<Tensor> ReferenceState::fock(const std::string& space) const {
72,376✔
419
  const auto itfound = m_fock.find(space);
72,376✔
420
  if (itfound != m_fock.end()) return itfound->second;
72,376✔
421

422
  std::unique_ptr<MoIndexTranslation> trans_ptr;
2,094✔
423
  try {
424
    trans_ptr.reset(new MoIndexTranslation(m_mo_ptr, space));
2,094✔
425
  } catch (const invalid_argument& inv) {
×
426
    throw invalid_argument("Invalid fock space string: " + std::string(inv.what()));
×
427
  }
×
428
  MoIndexTranslation& idxtrans = *trans_ptr;
2,094✔
429
  if (idxtrans.ndim() != 2) {
2,094✔
430
    throw invalid_argument("Invalid fock space string" + space +
×
431
                           ": Needs to consist of exactly 2 spaces.");
×
432
  }
433

434
  {
435
    RecordTime rec(m_timer, "import/fock/" + space);
2,094✔
436
    // Fock operator is symmetric wrt. index permutation and transforms
437
    // as the totally symmetric irrep
438
    const bool symmetric = true;
2,094✔
439
    std::shared_ptr<Symmetry> sym_ptr =
440
          make_symmetry_operator(m_mo_ptr, space, symmetric, "1");
4,188✔
441
    std::shared_ptr<Tensor> fock    = make_tensor_zero(sym_ptr);
2,094✔
442
    const HartreeFockSolution_i& hf = *m_hfsoln_ptr;
2,094✔
443

444
    auto fock_generator = [&idxtrans, &hf](
3,956✔
445
                                const std::vector<std::pair<size_t, size_t>>& block,
446
                                scalar_type* buffer) {
7,928✔
447
      // Strides, which should be used to write data to buffer
448
      // and size of the buffer.
449
      const std::vector<size_t> blength{block[0].second - block[0].first,
3,956✔
450
                                        block[1].second - block[1].first};
3,956✔
451
      const std::vector<size_t> strides{blength[1], 1};
3,956✔
452
      const size_t buffer_size = strides[0] * blength[0];
3,956✔
453

454
      std::vector<RangeMapping> mappings = idxtrans.map_range_to_hf_provider(block);
7,912✔
455
      for (RangeMapping& map : mappings) {
7,928✔
456
        const SimpleRange& bidx  = map.from();  // Indices of the MO subspace block
3,972✔
457
        const SimpleRange& hfidx = map.to();    // Indices translated to HF provider
3,972✔
458

459
        // Compute offset and remaining buffer size
460
        size_t offset = 0;
3,972✔
461
        for (size_t i = 0; i < 2; ++i) {
11,916✔
462
          offset += (bidx.axis(i).start() - block[i].first) * strides[i];
7,944✔
463
        }
464
        const size_t size = buffer_size - offset;
3,972✔
465
        hf.fock_ff(hfidx.axis(0).start(), hfidx.axis(0).end(),  //
3,972✔
466
                   hfidx.axis(1).start(), hfidx.axis(1).end(),  //
7,944✔
467
                   strides[0], strides[1], buffer + offset, size);
3,972✔
468
      }
469
    };
3,956✔
470
    fock->import_from(fock_generator, hf.conv_tol(), m_symmetry_check_on_import);
2,094✔
471
    fock->set_immutable();
2,094✔
472

473
    m_fock[space] = fock;
2,094✔
474
    return fock;
2,094✔
475
  }
2,094✔
476
}
2,094✔
477

478
std::shared_ptr<Tensor> ReferenceState::eri(const std::string& space) const {
198,608✔
479
  const auto itfound = m_eri.find(space);
198,608✔
480
  if (itfound != m_eri.end()) return itfound->second;
198,608✔
481

482
  std::unique_ptr<MoIndexTranslation> trans_ptr;
3,892✔
483
  try {
484
    trans_ptr.reset(new MoIndexTranslation(m_mo_ptr, space));
3,892✔
485
  } catch (const invalid_argument& inv) {
×
486
    throw invalid_argument("Invalid ERI space string: " + std::string(inv.what()));
×
487
  }
×
488
  if (trans_ptr->ndim() != 4) {
3,892✔
489
    throw invalid_argument("Invalid ERI space string" + space +
×
490
                           ": Needs to consist of exactly 4 spaces.");
×
491
  }
492

493
  // Branch off to three different ways to import the tensor, depending on the
494
  // capabilities of the HartreFockSolution_i and the block to be imported.
495
  RecordTime rec(m_timer, "import/eri/" + space);
3,892✔
496
  const std::vector<std::string>& ss = trans_ptr->subspaces();
3,892✔
497
  if (m_hfsoln_ptr->has_eri_phys_asym_ffff()) {
3,892✔
498
    m_eri[space] =
×
499
          import_eri_asym_direct(*m_hfsoln_ptr, *trans_ptr, m_symmetry_check_on_import);
×
500
  } else if (ss[2] == ss[3] || ss[0] == ss[1]) {
3,892✔
501
    m_eri[space] = import_eri_chem_then_asym_fast(*m_hfsoln_ptr, *trans_ptr,
5,420✔
502
                                                  m_symmetry_check_on_import);
5,420✔
503
  } else {
504
    m_eri[space] = import_eri_chem_then_asym(*m_hfsoln_ptr, *trans_ptr,
2,364✔
505
                                             m_symmetry_check_on_import);
2,364✔
506
  }
507
  m_eri[space]->set_immutable();
3,892✔
508
  return m_eri[space];
3,892✔
509
}
3,892✔
510

511
void ReferenceState::import_all() const {
22✔
512
  auto& ss = m_mo_ptr->subspaces;
22✔
513
  std::vector<std::string> ss_pairs;
22✔
514
  for (auto it1 = ss.begin(); it1 != ss.end(); ++it1) {
76✔
515
    for (auto it2 = it1; it2 != ss.end(); ++it2) {
150✔
516
      ss_pairs.push_back(*it1 + *it2);
96✔
517
      fock(*it1 + *it2);
96✔
518
    }
519
  }
520

521
  for (auto it1 = ss_pairs.begin(); it1 != ss_pairs.end(); ++it1) {
118✔
522
    for (auto it2 = it1; it2 != ss_pairs.end(); ++it2) {
378✔
523
      eri(*it1 + *it2);
282✔
524
    }
525
  }
526

527
  flush_hf_cache();
22✔
528
}
22✔
529

530
std::vector<std::string> ReferenceState::cached_fock_blocks() const {
×
531
  std::vector<std::string> ret;
×
532
  for (auto& kv : m_fock) ret.push_back(kv.first);
×
533
  return ret;
×
534
}
×
535

536
std::vector<std::string> ReferenceState::cached_eri_blocks() const {
×
537
  std::vector<std::string> ret;
×
538
  for (auto& kv : m_eri) ret.push_back(kv.first);
×
539
  return ret;
×
540
}
×
541

542
void ReferenceState::set_cached_fock_blocks(std::vector<std::string> newlist) {
×
543
  for (auto& item : newlist) fock(item);
×
544
  flush_hf_cache();
×
545

546
  std::sort(newlist.begin(), newlist.end());
×
547
  auto it = m_fock.begin();
×
548
  while (it != m_fock.end()) {
×
549
    if (std::binary_search(newlist.begin(), newlist.end(), it->first)) {
×
550
      ++it;
×
551
    } else {
552
      it = m_fock.erase(it);
×
553
    }
554
  }
555
}
×
556

557
void ReferenceState::set_cached_eri_blocks(std::vector<std::string> newlist) {
×
558
  for (auto& item : newlist) eri(item);
×
559
  flush_hf_cache();
×
560

561
  std::sort(newlist.begin(), newlist.end());
×
562
  auto it = m_eri.begin();
×
563
  while (it != m_eri.end()) {
×
564
    if (std::binary_search(newlist.begin(), newlist.end(), it->first)) {
×
565
      ++it;
×
566
    } else {
567
      it = m_eri.erase(it);
×
568
    }
569
  }
570
}
×
571

572
}  // namespace libadcc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc