MPQC  3.0.0-alpha
iters.h
1 
2 // cadf_iters.h
3 //
4 // Copyright (C) 2014 David Hollman
5 //
6 // Author: David Hollman
7 // Maintainer: DSH
8 // Created: Jan 10, 2014
9 //
10 // This file is part of the SC Toolkit.
11 //
12 // The SC Toolkit is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Library General Public License as published by
14 // the Free Software Foundation; either version 2, or (at your option)
15 // any later version.
16 //
17 // The SC Toolkit is distributed in the hope that it will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 // GNU Library General Public License for more details.
21 //
22 // You should have received a copy of the GNU Library General Public License
23 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
24 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
25 //
26 // The U.S. Government is granted a limited license as per AL 91-7.
27 //
28 
29 #ifndef _chemistry_qc_scf_cadf_iters_h
30 #define _chemistry_qc_scf_cadf_iters_h
31 
32 // standard library includes
33 #include <utility>
34 #include <type_traits>
35 #include <memory>
36 #include <iterator>
37 #include <atomic>
38 #include <unordered_set>
39 #include <unordered_map>
40 #include <limits>
41 
42 // boost includes
43 #include <boost/thread/thread.hpp>
44 #include <boost/functional/hash.hpp>
45 #include <boost/range.hpp>
46 #include <boost/range/irange.hpp>
47 #include <boost/range/join.hpp>
48 #include <boost/range/counting_range.hpp>
49 #include <boost/iterator/iterator_adaptor.hpp>
50 #include <boost/iterator/iterator_facade.hpp>
51 #include <boost/unordered_set.hpp>
52 
53 // MPQC includes
54 #include <util/misc/thread_timer.h>
55 #include <util/misc/iterators.h>
56 #include <chemistry/qc/scf/clhf.h>
57 
58 #include "iters_fwd.h"
59 #include "macros.h"
60 
61 #define DEFAULT_TARGET_BLOCK_SIZE 10000 // functions
62 #define DEFAULT_MAX_AM 7
63 
64 namespace sc {
65 
66 enum {
67  NotAssigned = -998,
68  NoLastIndex = -999,
69  NoMaximumBlockSize = -1001,
70  IndicesEnd = -1002
71 };
72 
73 using int_range = decltype(boost::counting_range(1, 2));
74 
75 struct ignored_argument { };
76 
77 template <typename T>
79  T* p_ = 0;
80 
81  public:
82 
83  OptionalRefParameter(T* p) : p_(p) { }
84 
85  OptionalRefParameter(const Ref<T>& rp) : p_(rp.nonnull() ? rp.pointer() : 0) { }
86 
87  operator T*() const { return p_; }
88  operator const T*() const { return p_; }
89 
90 };
91 
92 typedef enum {
93  NoRestrictions = 0,
94  SameCenter = 1,
95  SameAngularMomentum = 2,
96  Contiguous = 4
97 } BlockCompositionRequirement;
98 
99 //############################################################################//
100 
101 template <typename Iterator, typename DataContainer> struct BasisElementIteratorDereference;
102 
104 
105  public:
106 
108  int index,
109  GaussianBasisSet* basis,
110  GaussianBasisSet* dfbasis
111  ) : index(index),
112  basis(basis),
113  dfbasis(dfbasis)
114  { }
115 
116  BasisElementData() = delete;
117 
118  GaussianBasisSet* basis;
119  GaussianBasisSet* dfbasis;
120 
121  int index = NotAssigned;
122 
123  bool operator!=(const BasisElementData& other) const
124  {
125  return index != other.index;
126  }
127 
128  virtual ~BasisElementData() { }
129 
130  int center = NotAssigned;
131  int bfoff_in_atom = NotAssigned;
132 
133 };
134 
135 struct ShellData : public BasisElementData {
136 
137  template<typename Iterator> using with_iterator =
139 
140  ShellData(
141  int index,
142  GaussianBasisSet* basis,
143  GaussianBasisSet* dfbasis = 0,
144  int block_offset = NotAssigned
145  ) : BasisElementData(index, basis, dfbasis),
146  block_offset(NotAssigned)
147  {
148  init();
149  }
150 
151  template <typename Iterator>
152  ShellData(
153  const with_iterator<Iterator>& deref
154  ) : ShellData(deref.index, deref.basis, deref.dfbasis)
155  { }
156 
157  ShellData() : BasisElementData(NotAssigned, 0, 0) { }
158 
159 
160  int bfoff = NotAssigned;
161  int nbf = NotAssigned;
162  int atom_bfoff = NotAssigned;
163  int atom_shoff = NotAssigned;
164  int atom_nsh = NotAssigned;
165  int atom_nbf = NotAssigned;
166  int shoff_in_atom = NotAssigned;
167  int atom_last_function = NotAssigned;
168  int atom_last_shell = NotAssigned;
169  int last_function = NotAssigned;
170  int am = NotAssigned;
171  int ncontraction = NotAssigned;
172 
173  // Used when an auxiliary basis is set in the parent ShellIter. Otherwise, set to NotAssigned
174  // The atom_obs* versions make more sense in cases where the primary basis is the dfbasis
175  union { int atom_dfshoff = NotAssigned; int atom_obsshoff; };
176  union { int atom_dfbfoff = NotAssigned; int atom_obsbfoff; };
177  union { int atom_dfnbf = NotAssigned; int atom_obsnbf; };
178  union { int atom_dfnsh = NotAssigned; int atom_obsnsh; };
179  union { int atom_df_last_function = NotAssigned; int atom_obs_last_function; };
180  union { int atom_df_last_shell = NotAssigned; int atom_obs_last_shell; };
181 
182 
183  // Reserved for future use
184  int block_offset = NotAssigned;
185 
186  operator const int() const {
187  out_assert(basis, !=, 0);
188  out_assert(index, <, basis->nshell());
189  out_assert(index, >=, 0);
190  return index;
191  }
192 
193  static int max_index(GaussianBasisSet* basis) {
194  return basis->nshell();
195  }
196 
197  bool is_generally_contracted() {
198  return ncontraction > 1;
199  }
200 
201  protected:
202 
203  void init(){
204 
205 
206  if(not basis) {
207  index = NotAssigned;
208  return;
209  }
210  if(index == NotAssigned || index == basis->nshell()) return;
211  if(index >= basis->nshell() || index < 0) return;
212 
213  const auto& shell = basis->shell(index);
214  nbf = shell.nfunction();
215  am = shell.am(0);
216  ncontraction = shell.ncontraction();
217  bfoff = basis->shell_to_function(index);
218  center = basis->shell_to_center(index);
219  atom_shoff = basis->shell_on_center(center, 0);
220  atom_bfoff = basis->shell_to_function(atom_shoff);
221  atom_nsh = basis->nshell_on_center(center);
222  atom_nbf = basis->nbasis_on_center(center);
223  shoff_in_atom = index - atom_shoff;
224  bfoff_in_atom = bfoff - atom_bfoff;
225  atom_last_function = atom_bfoff + atom_nbf - 1;
226  last_function = bfoff + nbf - 1;
227  atom_last_shell = atom_shoff + atom_nsh - 1;
228  if(dfbasis != 0) {
229  atom_dfshoff = dfbasis->shell_on_center(center, 0);
230  atom_dfbfoff = dfbasis->shell_to_function(atom_dfshoff);
231  atom_dfnsh = dfbasis->nshell_on_center(center);
232  atom_dfnbf = dfbasis->nbasis_on_center(center);
233  atom_df_last_function = atom_dfbfoff + atom_dfnbf - 1;
234  atom_df_last_shell = atom_dfshoff + atom_dfnsh - 1;
235  }
236 
237  }
238 };
239 
241 
242  template<typename Iterator> using with_iterator =
244 
246  int idx,
247  GaussianBasisSet* basis,
248  GaussianBasisSet* dfbasis,
249  int block_offset = NotAssigned
250  ) : BasisElementData(idx, basis, dfbasis),
251  block_offset(block_offset)
252  {
253  init();
254  }
255 
256  template <typename Iterator>
258  const with_iterator<Iterator>& deref
259  ) : BasisFunctionData(deref.index, deref.basis, deref.dfbasis, deref.block_offset)
260  { }
261 
263  : BasisFunctionData(NotAssigned, 0, 0)
264  { }
265 
266  int shell_index = NotAssigned;
267  int shell_bfoff = NotAssigned;
268  union { int bfoff_in_shell = NotAssigned; int off; };
269  int atom_nbf = NotAssigned;
270  int atom_shoff = NotAssigned;
271  int block_offset = NotAssigned;
272  int atom_bfoff = NotAssigned;
273  int bfoff_in_block = NotAssigned;
274 
275  // Used when an auxiliary basis is set. Otherwise, set to NotAssigned.
276  // The atom_obs* versions make more sense in cases where the primary basis is the dfbasis
277  union { int atom_dfnbf = NotAssigned; int atom_obsnbf; };
278  union { int atom_dfshoff = NotAssigned; int atom_obsshoff; };
279  union { int atom_dfbfoff = NotAssigned; int atom_obsbfoff; };
280 
281  operator const int() const { return index; }
282 
283  static int max_index(GaussianBasisSet* basis) {
284  return basis->nbasis();
285  }
286 
287  protected:
288 
289  void init()
290  {
291  if(index == NotAssigned || index == basis->nbasis()) return;
292 
293  out_assert(basis, !=, 0);
294  out_assert(index, <, basis->nbasis());
295  out_assert(index, >=, 0);
296 
297  shell_index = basis->function_to_shell(index);
298  shell_bfoff = basis->shell_to_function(shell_index);
299  center = basis->shell_to_center(shell_index);
300  bfoff_in_shell = index - shell_bfoff;
301  atom_shoff = basis->shell_on_center(center, 0);
302  atom_bfoff = basis->shell_to_function(atom_shoff);
303  bfoff_in_atom = index - atom_bfoff;
304  atom_nbf = basis->nbasis_on_center(center);
305  if(dfbasis != 0){
306  atom_dfshoff = dfbasis->shell_on_center(center, 0);
307  atom_dfbfoff = dfbasis->shell_to_function(atom_dfshoff);
308  atom_dfnbf = dfbasis->nbasis_on_center(center);
309  }
310  if(block_offset != NotAssigned) {
311  bfoff_in_block = index - block_offset;
312  }
313  }
314 };
315 
316 template <typename DataContainer, typename Iterator>
317 struct BasisElementIteratorDereference : DataContainer {
318 
319  const Iterator iterator;
320 
321  template<typename... Args>
323  Iterator iter,
324  Args&&... args
325  ) : DataContainer(std::forward<Args>(args)...), iterator(iter)
326  { }
327 
328 };
329 
330 //############################################################################//
331 
332 template <typename DataContainer, typename Iterator=int_range::iterator>
334  : public boost::iterator_adaptor<
335  basis_element_iterator<DataContainer, Iterator>,
336  Iterator,
337  boost::use_default,
338  boost::bidirectional_traversal_tag,
339  BasisElementIteratorDereference<DataContainer, Iterator>
340  >
341 {
342 
343  public:
344 
346  typedef boost::iterator_adaptor<
347  self_t, Iterator, boost::use_default, boost::bidirectional_traversal_tag,
349  > super_t;
350  typedef typename DataContainer::template with_iterator<Iterator> reference;
351  typedef Iterator base_iterator;
352 
353 
354  protected:
355 
356  GaussianBasisSet* basis;
357  GaussianBasisSet* dfbasis;
358  int block_offset = NotAssigned;
359 
360  friend class boost::iterator_core_access;
361 
362  reference dereference() const {
363  auto base_spot = super_t::base();
364  return reference(
365  base_spot,
366  int(*base_spot),
367  basis, dfbasis,
368  block_offset
369  );
370  }
371 
372  public:
373 
374  basis_element_iterator() : basis(0), dfbasis(0) { }
375 
377  GaussianBasisSet* basis,
378  GaussianBasisSet* dfbasis,
379  Iterator iter,
380  int block_offset = NotAssigned
381  ) : basis_element_iterator::iterator_adaptor_(iter),
382  basis(basis),
383  dfbasis(dfbasis),
384  block_offset(block_offset)
385  { }
386 
387  operator Iterator() const {
388  return super_t::base();
389  }
390 
391 };
392 
393 template <typename DataContainer, typename Iterator=int_range::iterator>
395  : public boost::iterator_adaptor<
396  basis_element_with_value_iterator<DataContainer, Iterator>,
397  Iterator,
398  boost::use_default,
399  boost::bidirectional_traversal_tag,
400  BasisElementIteratorDereference<DataContainer, Iterator>
401  >
402 {
403 
404  public:
405 
407  typedef boost::iterator_adaptor<
408  self_t, Iterator, boost::use_default, boost::bidirectional_traversal_tag,
410  > super_t;
412 
413  basis_element_with_value_iterator() : basis(0), dfbasis(0) { }
414 
416  GaussianBasisSet* basis,
417  GaussianBasisSet* dfbasis,
418  Iterator iter,
419  int block_offset = NotAssigned
420  ) : super_t(iter),
421  basis(basis),
422  dfbasis(dfbasis),
423  block_offset(block_offset)
424  { }
425 
427  {
429  return result_type(basis, dfbasis, super_t::base(), block_offset);
430  }
431 
432  private:
433 
434  GaussianBasisSet* basis;
435  GaussianBasisSet* dfbasis;
436  int block_offset = NotAssigned;
437 
438  friend class boost::iterator_core_access;
439 
440  reference dereference() const {
441  auto base_spot = super_t::base();
442  return reference(
443  base_spot,
444  (*base_spot).index,
445  (*base_spot).value,
446  basis, dfbasis,
447  block_offset
448  );
449  }
450 };
451 
452 template<typename Iterator=int_range::iterator>
454 
455 template<typename Iterator=int_range::iterator>
457 
458 
459 //############################################################################//
460 
461 template<typename DataContainer, typename Iterator=int_range::iterator>
462  using range_of = decltype(boost::make_iterator_range(
465  ));
466 
467 template<typename DataContainer, typename Iterator=int_range::iterator>
468  using valued_range_of = decltype(boost::make_iterator_range(
471  ));
472 
473 template <typename DataContainer>
474 inline range_of<DataContainer>
475 basis_element_range(
476  GaussianBasisSet* basis,
478  int first_index,
479  int last_index,
480  int block_offset = NotAssigned
481 )
482 {
483  const int actual_last_index = last_index == NoLastIndex ?
484  DataContainer::max_index(basis) - 1 : last_index;
485  auto base_range = boost::counting_range(first_index, actual_last_index + 1);
486  return boost::make_iterator_range(
487  basis_element_iterator<DataContainer>(basis, dfbasis, base_range.begin(), block_offset),
488  basis_element_iterator<DataContainer>(basis, dfbasis, base_range.end(), block_offset)
489  );
490 
491 }
492 
493 template <typename DataContainer>
494 inline range_of<DataContainer>
495 basis_element_range(
496  GaussianBasisSet* basis,
497  int last_index
498 )
499 {
500  return basis_element_range<DataContainer>(basis, 0, 0, last_index);
501 }
502 
503 template <typename DataContainer>
504 inline range_of<DataContainer>
505 basis_element_range(
506  GaussianBasisSet* basis,
507  GaussianBasisSet* dfbasis = 0,
508  int last_index = NoLastIndex
509 )
510 {
511  return basis_element_range<DataContainer>(basis, dfbasis, 0, last_index);
512 }
513 
514 template <typename DataContainer>
515 inline range_of<DataContainer>
516 basis_element_range(
517  GaussianBasisSet* basis,
518  int first_index,
519  int last_index
520 )
521 {
522  return basis_element_range<DataContainer>(basis, 0, first_index, last_index);
523 }
524 
525 template <typename Iterable>
526 typename boost::enable_if<
527  boost::is_convertible<typename Iterable::value_type, int>,
528  range_of<ShellData, typename Iterable::iterator>
529 >::type
531  const Iterable& index_iterable,
532  GaussianBasisSet* basis,
533  GaussianBasisSet* dfbasis=0
534 )
535 {
536  return boost::make_iterator_range(
537  basis_element_iterator<ShellData, typename Iterable::iterator>(basis, dfbasis, index_iterable.begin()),
538  basis_element_iterator<ShellData, typename Iterable::iterator>(basis, dfbasis, index_iterable.end())
539  );
540 }
541 
542 template <typename Iterable>
543 typename boost::enable_if<
544  boost::is_convertible<typename Iterable::value_type, int>,
545  range_of<BasisFunctionData, typename Iterable::iterator const>
546 >::type
548  const Iterable& index_iterable,
549  GaussianBasisSet* basis,
550  GaussianBasisSet* dfbasis=0
551 )
552 {
553  return boost::make_iterator_range(
554  basis_element_iterator<BasisFunctionData, typename Iterable::iterator const>(basis, dfbasis, index_iterable.begin()),
555  basis_element_iterator<BasisFunctionData, typename Iterable::iterator const>(basis, dfbasis, index_iterable.end())
556  );
557 }
558 
559 template <typename Iterator>
560 inline range_of<ShellData, Iterator>
562  const ShellData::with_iterator<Iterator>& begin,
563  const ShellData::with_iterator<Iterator>& end
564 )
565 {
566  return boost::make_iterator_range(
567  basis_element_iterator<ShellData, Iterator>(begin.basis, begin.dfbasis, begin.iterator),
568  basis_element_iterator<ShellData, Iterator>(begin.basis, begin.dfbasis, end.iterator)
569  );
570 }
571 
572 template <typename Iterator, typename Iterator2>
573 inline range_of<ShellData, Iterator>
575  const ShellData::with_iterator<Iterator>& begin,
576  const Iterator2& end
577 )
578 {
579  return boost::make_iterator_range(
580  basis_element_iterator<ShellData, Iterator>(begin.basis, begin.dfbasis, begin.iterator),
581  basis_element_iterator<ShellData, Iterator>(begin.basis, begin.dfbasis, Iterator(end))
582  );
583 }
584 
585 template <typename Iterator>
586 inline range_of<ShellData, Iterator>
588  const Iterator& begin,
589  const Iterator& end,
590  GaussianBasisSet* basis,
591  GaussianBasisSet* dfbasis
592 )
593 {
594  return boost::make_iterator_range(
595  basis_element_iterator<ShellData, Iterator>(basis, dfbasis, begin),
596  basis_element_iterator<ShellData, Iterator>(basis, dfbasis, end)
597  );
598 }
599 
600 template <typename Iterator>
601 inline range_of<BasisFunctionData, Iterator>
603  const BasisFunctionData::with_iterator<Iterator>& begin,
604  const BasisFunctionData::with_iterator<Iterator>& end
605 )
606 {
607  return boost::make_iterator_range(
608  basis_element_iterator<BasisFunctionData, Iterator>(begin.basis, begin.dfbasis, begin.iterator),
609  basis_element_iterator<BasisFunctionData, Iterator>(begin.basis, begin.dfbasis, end.iterator)
610  );
611 }
612 
613 template <typename... Args>
614 inline range_of<ShellData>
615 shell_range(GaussianBasisSet* basis, Args&&... args) {
616  return basis_element_range<ShellData>(basis, std::forward<Args>(args)...);
617 }
618 
619 template <typename... Args>
620 inline range_of<BasisFunctionData>
621 function_range(GaussianBasisSet* basis, Args&&... args) {
622  return basis_element_range<BasisFunctionData>(basis, std::forward<Args>(args)...);
623 }
624 
625 inline range_of<BasisFunctionData>
626 function_range(const ShellData& ish) {
627  return basis_element_range<BasisFunctionData>(
628  ish.basis, ish.dfbasis,
629  ish.bfoff, ish.last_function
630  );
631 }
632 
633 //############################################################################//
634 
635 template<typename Iterator> class ShellBlockSkeleton;
636 
637 template<typename Range = range_of<ShellData>>
639 
640  public:
641 
642  typedef Range shell_range_t;
643 
644  protected:
645 
646  void init();
647 
648  public:
649 
650  bool contiguous_ = false;
651  Range shell_range;
652  ShellData first_shell;
653  ShellData last_shell;
654  GaussianBasisSet* basis;
655  GaussianBasisSet* dfbasis;
656  int restrictions;
657 
658  ShellBlockData() : restrictions(NotAssigned) { }
659 
661 
662  explicit ShellBlockData(
663  GaussianBasisSet* basis,
664  GaussianBasisSet* dfbasis = 0
665  ) : ShellBlockData(
666  sc::shell_range(basis, dfbasis, 0, basis->nshell() - 1),
667  basis->nshell(), basis->nbasis(),
668  NoRestrictions
669  )
670  { contiguous_ = true; }
671 
672  // Construct a one-shell shell block
673  ShellBlockData(const ShellData& ish)
674  : ShellBlockData(
675  sc::shell_range(ish.basis, ish.dfbasis, ish, ish),
676  1, ish.nbf, SameCenter|SameAngularMomentum
677  )
678  { }
679 
681  Range sh_range,
682  int nshell, int nbf,
683  int requirements
684  ) : shell_range(sh_range),
685  first_shell(*shell_range.begin()),
686  last_shell(shell_range.begin() == shell_range.end() ?
687  *shell_range.end() : *(--shell_range.end())
688  ),
689  basis(first_shell.basis),
690  dfbasis(first_shell.dfbasis),
691  nshell(nshell), nbf(nbf),
692  restrictions(requirements)
693  { init(); }
694 
695  static ShellBlockData atom_block(int atom,
696  GaussianBasisSet* basis,
697  GaussianBasisSet* dfbasis = 0
698  )
699  {
700  const int atom_shoff = basis->shell_on_center(atom, 0);
701  const int atom_nsh = basis->nshell_on_center(atom);
702  const int atom_nbf = basis->nbasis_on_center(atom);
703  return ShellBlockData<Range>(
704  sc::shell_range(basis, dfbasis, atom_shoff, atom_shoff + atom_nsh-1),
705  atom_nsh, atom_nbf, SameCenter
706  );
707  }
708 
709  template <typename OtherRange>
710  auto operator+(const ShellBlockData<OtherRange>& other)
711  -> ShellBlockData<decltype(boost::join(this->shell_range, other.shell_range))>
712  {
713  typedef ShellBlockData<decltype(boost::join(this->shell_range, other.shell_range))> result_t;
714  if(basis != other.basis) {
715  throw ProgrammingError("Can't combine blocks with different basis sets", __FILE__, __LINE__);
716  }
717  if(dfbasis != other.dfbasis) {
718  throw ProgrammingError("Can't combine blocks with different auxiliary basis sets", __FILE__, __LINE__);
719  }
720  int new_reqs = NoRestrictions;
721  if(
722  (restrictions & SameCenter) and (other.restrictions & SameCenter)
723  and this->center == other.center
724  ) {
725  new_reqs |= SameCenter;
726  }
727  if(
728  (restrictions & SameAngularMomentum) and (other.restrictions & SameAngularMomentum)
729  and basis->shell(first_shell).am() == basis->shell(first_shell).am()
730  ) {
731  new_reqs |= SameAngularMomentum;
732  }
733  if(
734  (restrictions & Contiguous) and (other.restrictions & Contiguous)
735  and (int)last_shell + 1 == (int)other.first_shell
736  ) {
737  new_reqs |= Contiguous;
738  }
739 
740  return result_t(
741  boost::join(this->shell_range, other.shell_range),
742  nshell + other.nshell, nbf + other.nbf,
743  new_reqs
744  );
745 
746  }
747 
748  int nbf;
749  int bfoff;
750  int nshell;
751  int last_function;
752  int center = NotAssigned;
753  int atom_bfoff = NotAssigned;
754  int atom_shoff = NotAssigned;
755  int atom_nsh = NotAssigned;
756  int atom_nbf = NotAssigned;
757  int bfoff_in_atom = NotAssigned;
758  int shoff_in_atom = NotAssigned;
759  int atom_last_function = NotAssigned;
760  int atom_last_shell = NotAssigned;
761 
762  // Used when an auxiliary basis is set in the parent ShellIter. Otherwise, set to NotAssigned
763  // The atom_obs* versions make more sense in cases where the primary basis is the dfbasis
764  union { int atom_dfshoff = NotAssigned; int atom_obsshoff; };
765  union { int atom_dfbfoff = NotAssigned; int atom_obsbfoff; };
766  union { int atom_dfnbf = NotAssigned; int atom_obsnbf; };
767  union { int atom_dfnsh = NotAssigned; int atom_obsdfnsh; };
768  union { int atom_df_last_function = NotAssigned; int atom_obs_last_function; };
769  union { int atom_df_last_shell = NotAssigned; int atom_obs_last_shell; };
770 
771  static int max_index(const GaussianBasisSet* basis) { return basis->nshell(); }
772 
773  bool is_contiguous() const { return contiguous_; }
774 };
775 
776 template<typename Range = range_of<ShellData>>
777 class ShellBlockSkeleton {
778 
779  private:
780 
781  int nshell;
782  int nbf;
783 
784 
785  friend class ShellBlockData<Range>;
786 
787  public:
788 
789  Range shell_range;
790  int first_index;
791  int restrictions;
792 
793  ShellBlockSkeleton() : nshell(0), nbf(0), first_index(NotAssigned) { }
794 
795  static ShellBlockSkeleton<Range> end_skeleton() { return ShellBlockSkeleton(); }
796 
797  ShellBlockSkeleton(const ShellBlockData<Range>& shbd)
798  : shell_range(shbd.shell_range), nshell(shbd.nshell), nbf(shbd.nbf),
799  restrictions(shbd.restrictions),
800  first_index((*shell_range.begin()).index)
801  { }
802 
803  ShellBlockSkeleton(Range range, int nsh, int nbas, int restrictions)
804  : shell_range(range), nshell(nsh), nbf(nbas),
805  restrictions(restrictions),
806  first_index((*shell_range.begin()).index)
807  { }
808 
809 
810  template<typename OtherRange>
811  bool operator==(const ShellBlockSkeleton<OtherRange>& other) const
812  {
813  return !(this->operator!=(other));
814  }
815 
816  template<typename OtherRange>
817  bool operator!=(const ShellBlockSkeleton<OtherRange>& other) const
818  {
819  return first_index != other.first_index
820  or nshell != other.nshell
821  or nbf != other.nbf;
822  }
823 
824  template<typename OtherRange>
825  bool operator<(const ShellBlockSkeleton<OtherRange>& other) const
826  {
827  const int my_index = (*shell_range.begin()).index;
828  const int other_index = (*other.shell_range.begin()).index;
829  if(my_index < other_index) return true;
830  else if(my_index > other_index) return false;
831  else if(nshell < other.nshell) return true;
832  else if(nshell > other.nshell) return false;
833  else if(nbf < other.nbf) return true;
834  else return false;
835  }
836 
837 };
838 
839 template<typename ShellIterator, typename ShellRange=range_of<ShellData, ShellIterator>>
841  : public boost::iterator_facade<
842  shell_block_iterator<ShellIterator, ShellRange>,
843  ShellBlockSkeleton<range_of<ShellData, ShellIterator>>,
844  boost::forward_traversal_tag,
845  ShellBlockData<range_of<ShellData, ShellIterator>>
846  >
847 {
848  public:
849 
852 
853  private:
854 
855  int target_size;
856  int restrictions;
857 
858  ShellRange all_shells;
859  ShellBlockSkeleton<ShellRange> current_skeleton;
860 
861  void init();
862  void init_from_spot(const decltype(all_shells.begin())& start_spot);
863 
864  bool contiguous_ = false;
865 
866  friend class boost::iterator_core_access;
867 
868  value_reference dereference() const {
869  auto rv = value_reference(
870  current_skeleton
871  );
872  rv.contiguous_ = contiguous_;
873  return rv;
874  }
875 
876  template <typename OtherIterator>
877  bool equal(shell_block_iterator<OtherIterator> const& other) const
878  {
879  return current_skeleton == other.current_skeleton;
880  }
881 
882  void increment()
883  {
884  // Don't allow incrementing of end iterator
885  assert(current_skeleton.first_index != NotAssigned);
886  auto&& new_begin = current_skeleton.shell_range.end();
887  auto&& new_end = all_shells.end();
888  all_shells = boost::make_iterator_range(new_begin, new_end);
889  init();
890  }
891 
892  public:
893 
894  GaussianBasisSet* basis;
895  GaussianBasisSet* dfbasis;
896 
897  shell_block_iterator() : basis(0), dfbasis(0), restrictions(0), target_size(0) { }
898 
900  const ShellRange& all_shells_in,
901  int requirements = SameCenter,
902  int target_size = DEFAULT_TARGET_BLOCK_SIZE
903  ) : basis(all_shells_in.begin()->basis),
904  dfbasis(all_shells_in.begin()->dfbasis),
905  restrictions(requirements),
906  target_size(target_size),
907  all_shells(all_shells_in)
908  {
909  init();
910  }
911 
913  const ShellIterator& index_begin,
914  const ShellIterator& index_end,
915  GaussianBasisSet* basis,
916  GaussianBasisSet* dfbasis = 0,
917  int requirements = SameCenter,
918  int target_size = DEFAULT_TARGET_BLOCK_SIZE
919  ) : basis(basis),
920  dfbasis(dfbasis),
921  restrictions(requirements),
922  target_size(target_size),
923  all_shells(shell_range(index_begin, index_end, basis, dfbasis))
924  {
925  init();
926  }
927 
929  GaussianBasisSet* basis,
930  GaussianBasisSet* dfbasis = 0,
931  int first_index = 0,
932  int last_index = NoLastIndex,
933  int requirements = SameCenter,
934  int target_size = DEFAULT_TARGET_BLOCK_SIZE
935  ) : basis(basis),
936  dfbasis(dfbasis),
937  restrictions(requirements),
938  target_size(target_size),
939  all_shells(shell_range(
940  basis, dfbasis, first_index,
941  last_index == NoLastIndex ? ShellData::max_index(basis) - 1 : last_index
942  ))
943  {
944  init();
945  }
946 
947  static shell_block_iterator end_of_basis(GaussianBasisSet* basis) {
948  return shell_block_iterator(basis, 0, basis->nshell(), basis->nshell()-1);
949  }
950 
951  static shell_block_iterator end_with_last_index(int index, GaussianBasisSet* basis) {
952  // The end iterator needs to start after the "last_index" that is passed in. But
953  // because last_index gets incremented before being made into a counting range,
954  // the last_index parameter of the end iterator should be one less than the first.
955  // This is counterintuitive, but not meant to be exposed to the outside world anyway.
956  return shell_block_iterator(basis, 0, index + 1, index);
957  }
958 
959  static shell_block_iterator end_of_range(
960  const ShellRange& all_shells_in,
961  int requirements = SameCenter,
962  int target_size = DEFAULT_TARGET_BLOCK_SIZE
963  )
964  {
965  return shell_block_iterator(
966  boost::make_iterator_range(
967  all_shells_in.end(),
968  all_shells_in.end()
969  ),
970  requirements, target_size
971  );
972  }
973 
974  const self_type& requiring(int restr) {
975  restrictions = restr;
976  init();
977  return *this;
978  }
979 
980 };
981 
982 // Type alias specialization for shell blocks
983 template<typename Iterator=int_range::iterator, typename Range=range_of<ShellData, Iterator>>
984  using range_of_shell_blocks = decltype(boost::make_iterator_range(
987  ));
988 
989 inline range_of_shell_blocks<>
991  GaussianBasisSet* basis,
993  int first_index = 0,
994  int last_index = NoLastIndex,
995  int reqs = SameCenter,
996  int target_size = DEFAULT_TARGET_BLOCK_SIZE
997 )
998 {
1000  const int actual_last_index = last_index == NoLastIndex ?
1001  basis->nshell() - 1 : last_index;
1002  return boost::make_iterator_range(
1003  result_t(basis, dfbasis, first_index,
1004  actual_last_index, reqs, target_size
1005  ),
1006  result_t::end_with_last_index(actual_last_index, basis)
1007  );
1008 
1009 }
1010 
1011 inline range_of_shell_blocks<>
1012 shells_blocked_by_atoms(
1013  GaussianBasisSet* basis,
1014  OptionalRefParameter<GaussianBasisSet> dfbasis = 0
1015 )
1016 {
1017  return shell_block_range(basis, dfbasis, 0, NoLastIndex, SameCenter, NoMaximumBlockSize);
1018 }
1019 
1020 inline std::ostream&
1021 operator << (std::ostream& out, const ShellData& ish)
1022 {
1023  out << "ShellData: {"
1024  << "index = " << ish.index << ", ";
1025  out << "center = ";
1026  if(ish.center == NotAssigned) out << "NotAssigned, ";
1027  else out << ish.center << ", ";
1028  out << "bfoff = ";
1029  if(ish.bfoff == NotAssigned) out << "NotAssigned, ";
1030  else out << ish.bfoff << ", ";
1031  out << "nbf = ";
1032  if(ish.nbf == NotAssigned) out << "NotAssigned, ";
1033  else out << ish.nbf << ", ";
1034  out << "last_function = ";
1035  if(ish.nbf == NotAssigned) out << "NotAssigned";
1036  else out << ish.last_function;
1037  out << "}";
1038  return out;
1039 }
1040 
1041 template<typename Range>
1042 std::ostream&
1043 operator << (std::ostream& out, const ShellBlockData<Range>& blk)
1044 {
1045  auto write_if_assigned = [&](int val) {
1046  if(val == NotAssigned) out << "NotAssigned";
1047  else out << val;
1048  };
1049  out << "\nShellBlockData: {"
1050  << "\n nshell: " << blk.nshell
1051  << ", nbf: " << blk.nbf;
1052  out << "\n shoff: ";
1053  write_if_assigned(blk.first_shell.index);
1054  out << ", bfoff: ";
1055  write_if_assigned(blk.bfoff);
1056  out << "\n (basis nshell = "
1057  << blk.basis->nshell()
1058  << ", nbf = "
1059  << blk.basis->nbasis()
1060  << ")"
1061  << "\n center: ";
1062  write_if_assigned(blk.center);
1063  out << "\n shells:";
1064  for(auto sh : shell_range(blk)) {
1065  out << "\n " << sh;
1066  }
1067  out << "\n}" << std::endl;
1068  return out;
1069 }
1070 
1071 //############################################################################//
1072 
1073 inline range_of<BasisFunctionData>
1074 iter_functions_on_center(
1075  const Ref<GaussianBasisSet>& basis,
1076  int center,
1077  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
1078 )
1079 {
1080  const int shoff = basis->shell_on_center(center, 0);
1081  const int bfoff = basis->shell_to_function(shoff);
1082  return function_range(
1083  basis, dfbasis,
1084  bfoff, bfoff + basis->nbasis_on_center(center) - 1
1085  );
1086 }
1087 
1088 inline range_of<ShellData>
1089 iter_shells_on_center(
1090  GaussianBasisSet* basis,
1091  int center,
1092  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
1093 )
1094 {
1095  const int shoff = basis->shell_on_center(center, 0);
1096  return shell_range(
1097  basis, dfbasis,
1098  shoff, shoff + basis->nshell_on_center(center) - 1
1099  );
1100 }
1101 
1102 template<typename DataContainer, typename Iterator=int_range::iterator>
1103 using joined_range_of = decltype(
1104  boost::join(
1105  range_of<DataContainer, Iterator>(),
1106  range_of<DataContainer, Iterator>()
1107  )
1108 );
1109 
1110 inline joined_range_of<ShellData>
1111 iter_shells_on_centers(
1112  const Ref<GaussianBasisSet>& basis,
1113  int center1,
1114  int center2,
1115  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
1116 )
1117 {
1118  const int shoff1 = basis->shell_on_center(center1, 0);
1119  const int shoff2 = basis->shell_on_center(center2, 0);
1120  if(center1 == center2) {
1121  return boost::join(
1122  shell_range(
1123  basis, dfbasis,
1124  shoff1, shoff1 + basis->nshell_on_center(center1) - 1
1125  ),
1126  // An empty range
1127  shell_range(
1128  basis, dfbasis,
1129  shoff1, shoff1 - 1
1130  )
1131  );
1132  }
1133  else {
1134  return boost::join(
1135  shell_range(
1136  basis, dfbasis,
1137  shoff1, shoff1 + basis->nshell_on_center(center1) - 1
1138  ),
1139  shell_range(
1140  basis, dfbasis,
1141  shoff2, shoff2 + basis->nshell_on_center(center2) - 1
1142  )
1143  );
1144  }
1145 }
1146 
1147 inline range_of_shell_blocks<>
1148 iter_shell_blocks_on_center(
1149  const Ref<GaussianBasisSet>& basis,
1150  int center,
1151  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0,
1152  int reqs = SameCenter,
1153  int target_size = DEFAULT_TARGET_BLOCK_SIZE
1154 )
1155 {
1156  const int shoff = basis->shell_on_center(center, 0);
1157  return shell_block_range(
1158  basis, dfbasis,
1159  shoff,
1160  shoff + basis->nshell_on_center(center) - 1,
1161  reqs|SameCenter,
1162  target_size
1163  );
1164 }
1165 
1166 //############################################################################//
1167 
1169 
1170  template<typename Iterator> using with_iterator =
1172 
1173  double value;
1174 
1176  int index,
1177  double value,
1178  GaussianBasisSet* basis,
1179  GaussianBasisSet* dfbasis,
1180  int block_offset = NotAssigned
1181  ) : ShellData(index, basis, dfbasis), value(value)
1182  { }
1183 
1185  const ShellData& ish,
1186  double value
1187  ) : ShellData(ish), value(value)
1188  { }
1189 
1190 };
1191 
1192 template <typename Iterator>
1193 inline valued_range_of<ShellDataWithValue, Iterator>
1194 shell_range(
1197 )
1198 {
1199  return boost::make_iterator_range(
1200  basis_element_with_value_iterator<ShellDataWithValue, Iterator>(begin.basis, begin.dfbasis, begin.iterator),
1201  basis_element_with_value_iterator<ShellDataWithValue, Iterator>(begin.basis, begin.dfbasis, end.iterator)
1202  );
1203 }
1204 
1205 template <typename Iterator, typename Iterator2>
1206 inline valued_range_of<ShellData, Iterator>
1207 shell_range(
1208  const ShellDataWithValue::with_iterator<Iterator>& begin,
1209  const Iterator2& end
1210 )
1211 {
1212  return boost::make_iterator_range(
1213  basis_element_with_value_iterator<ShellDataWithValue, Iterator>(begin.basis, begin.dfbasis, begin.iterator),
1214  basis_element_with_value_iterator<ShellDataWithValue, Iterator>(begin.basis, begin.dfbasis, Iterator(end))
1215  );
1216 }
1217 
1218 template <typename Iterator>
1219 inline valued_range_of<ShellData, Iterator>
1220 shell_range(
1221  const basis_element_with_value_iterator<ShellData, Iterator>& begin,
1222  const basis_element_with_value_iterator<ShellData, Iterator>& end
1223 )
1224 {
1225  return boost::make_iterator_range(begin, end);
1226 }
1227 
1230 
1231  int index;
1232  mutable double value;
1233 
1234  ShellIndexWithValue() = default;
1235 
1236  ShellIndexWithValue(int index)
1237  : index(index), value(0.0)
1238  { }
1239 
1240  ShellIndexWithValue(int index, double value)
1241  : index(index), value(value)
1242  { }
1243 
1244  bool operator==(const ShellIndexWithValue& other) {
1245  return index == other.index;
1246  }
1247 
1248  operator int() const {
1249  return index;
1250  }
1251 };
1252 
1253 namespace detail {
1254  template<typename T>
1255  struct hash_;
1256 
1257  template<>
1259 
1260  std::size_t operator()(const ShellIndexWithValue& v) const {
1261  return std::hash<int>()(v.index);
1262  }
1263 
1264  };
1265 
1266  struct index_equal_ {
1267  bool operator()(const ShellIndexWithValue& a, const ShellIndexWithValue& b) const {
1268  return a.index == b.index;
1269  }
1270  };
1271 }
1272 
1273 
1274 template<typename Iterable>
1275 range_of_shell_blocks<typename Iterable::const_iterator>
1277  const Iterable& shlist, GaussianBasisSet* basis, GaussianBasisSet* dfbasis = 0,
1278  int requirements=SameCenter,
1279  int target_size=DEFAULT_TARGET_BLOCK_SIZE
1280 )
1281 {
1282  return boost::make_iterator_range(
1284  shlist.begin(), shlist.end(),
1285  basis, dfbasis, requirements, target_size
1286  ),
1288  shlist.end(), shlist.end(),
1289  basis, dfbasis, requirements, target_size
1290  )
1291  );
1292 }
1293 
1294 template <typename ShellRange>
1295 inline range_of_shell_blocks<typename ShellRange::iterator::base_iterator, ShellRange>
1297  const ShellBlockData<ShellRange>& iblk,
1298  int requirements=SameCenter,
1299  int target_size=DEFAULT_TARGET_BLOCK_SIZE
1300 )
1301 {
1302  return boost::make_iterator_range(
1303  shell_block_iterator<typename ShellRange::iterator::base_iterator, ShellRange>(
1304  iblk.shell_range, requirements|iblk.restrictions, target_size
1305  ),
1306  shell_block_iterator<typename ShellRange::iterator::base_iterator, ShellRange>::end_of_range(
1307  iblk.shell_range, requirements|iblk.restrictions, target_size
1308  )
1309  );
1310 }
1311 
1312 //############################################################################//
1313 
1315 
1317  std::unordered_map<int, std::vector<skeleton_t>> lists_for_restrictions_;
1318 
1319  public:
1320 
1322 
1323  // NOTE: idxlist must be sorted!
1324  template<typename Iterable>
1326  const Iterable& idxlist,
1327  GaussianBasisSet* basis,
1328  GaussianBasisSet* dfbasis,
1329  std::set<int> restriction_sets
1330  )
1331  {
1332  for(auto restrictions : restriction_sets) {
1333  lists_for_restrictions_.emplace(
1334  std::piecewise_construct,
1335  std::forward_as_tuple(restrictions),
1336  std::forward_as_tuple(0)
1337  );
1338  for(auto&& iblk : shell_block_range(idxlist, basis, dfbasis, Contiguous|restrictions, NoMaximumBlockSize)) {
1339  lists_for_restrictions_[restrictions].emplace_back(
1340  shell_range(basis, dfbasis, (int)iblk.first_shell, (int)iblk.last_shell),
1341  iblk.nshell, iblk.nbf, restrictions
1342  );
1343  }
1344  }
1345 
1346  }
1347 
1348  const std::vector<ShellBlockSkeleton<>>&
1349  list_with_restrictions(int restrictions) const {
1350  assert(lists_for_restrictions_.find(restrictions) != lists_for_restrictions_.end());
1351  return lists_for_restrictions_.at(restrictions);
1352  }
1353 
1354  const std::vector<ShellBlockSkeleton<>>&
1355  operator[](int restrictions) const {
1356  return list_with_restrictions(restrictions);
1357  }
1358 
1359 };
1360 
1361 
1362 
1363 } // end namespace sc
1364 
1365 #include "cadf_iters_impl.h"
1366 
1367 #endif /* _chemistry_qc_scf_cadf_iters_h */
sc::BasisElementIteratorDereference
Definition: iters.h:101
sc::OptionalRefParameter
Definition: iters.h:78
ShellBlockIterator
Definition: cadf_attic.h:479
sc::ShellBlockSkeleton
Definition: iters.h:635
sc::GaussianBasisSet::shell_to_center
int shell_to_center(int ishell) const
Return the center on which the given shell is located.
Definition: gaussbas.h:512
sc::GaussianBasisSet::nshell
unsigned int nshell() const
Return the number of shells.
Definition: gaussbas.h:500
sc::Ref::nonnull
bool nonnull() const
Return !null().
Definition: ref.h:641
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::basis_element_with_value_iterator
Definition: iters.h:394
sc::detail::index_equal_
Definition: iters.h:1266
sc::ShellData
Definition: iters.h:135
shell_block_iterable
Definition: cadf_attic.h:406
sc::ignored_argument
Definition: iters.h:75
sc::Ref::pointer
T * pointer() const
Returns a pointer the reference counted object.
Definition: ref.h:413
sc::GaussianBasisSet::nbasis_on_center
unsigned int nbasis_on_center(int icenter) const
Return the number of basis functions on the given center.
sc::ShellBlockData
Definition: iters.h:638
sc::ShellDataWithValue
Definition: iters.h:1168
sc::ContiguousShellBlockList
Definition: iters.h:1314
sc::GaussianBasisSet::nbasis
unsigned int nbasis() const
Return the number of basis functions.
Definition: gaussbas.h:514
ShellData
Definition: cadf_attic.h:100
sc::GaussianBasisSet::function_to_shell
int function_to_shell(int i) const
Return the shell to which the given function belongs.
sc::operator<<
std::vector< unsigned int > operator<<(const GaussianBasisSet &B, const GaussianBasisSet &A)
computes a map from basis functions in A to the equivalent basis functions in B.
sc::detail::hash_
Definition: iters.h:1255
ShellBlockIterator::operator!=
bool operator!=(const ShellBlockIterator &other) const
Note: not a full != operator; just for the purposes of range iteration! Two blocks with the same firs...
Definition: cadf_attic.h:549
sc::other
SpinCase1 other(SpinCase1 S)
given 1-spin return the other 1-spin
shell_iterable
Definition: cadf_attic.h:330
sc::GaussianBasisSet::shell_on_center
int shell_on_center(int icenter, int shell) const
Return the overall shell number, given a center and the shell number on that center.
sc::GaussianShell::am
const std::vector< unsigned int > & am() const
The angular momenta of contractions.
Definition: gaussshell.h:199
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition: scexception.h:92
mpqc::ci::operator<
bool operator<(const Space< Spin > &a, const Space< Spin > &b)
Compare two spaces by their rank.
Definition: subspace.hpp:44
sc::shell_block_iterator
Definition: iters.h:840
sc::GaussianBasisSet
The GaussianBasisSet class is used describe a basis set composed of atomic gaussian orbitals.
Definition: gaussbas.h:141
sc::BasisElementData
Definition: iters.h:103
sc::basis_element_iterator
Definition: iters.h:333
function_iterable
Definition: cadf_attic.h:375
sc::GaussianBasisSet::shell_to_function
int shell_to_function(int i) const
Return the number of the first function in the given shell.
Definition: gaussbas.h:540
sc::BasisFunctionData
Definition: iters.h:240
sc::operator==
bool operator==(const Atom &a, const Atom &b)
sc::GaussianBasisSet::nshell_on_center
int nshell_on_center(int icenter) const
Return the number of shells on the given center.
sc::GaussianBasisSet::shell
const Shell & shell(int i) const
Return a reference to Shell number i.
Definition: gaussbas.h:558
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::ShellIndexWithValue
binds an integer index + real annotation, e.g. Shell index + associated operator norm
Definition: iters.h:1229

Generated at Sun Jan 26 2020 23:23:59 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.