MPQC  3.0.0-alpha
cadf_attic.h
1 //
2 // cadf_attic.h
3 //
4 // Copyright (C) 2014 David Hollman
5 //
6 // Author: David Hollman
7 // Maintainer: DSH
8 // Created: Jan 13, 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 CADF_ATTIC_H_
30 #define CADF_ATTIC_H_
31 
32 #include <boost/preprocessor/seq.hpp>
33 #include <boost/preprocessor/logical/not.hpp>
34 #include <boost/preprocessor/control/if.hpp>
35 #include <boost/preprocessor/tuple/rem.hpp>
36 #include <boost/preprocessor/facilities/empty.hpp>
37 #include <boost/preprocessor/punctuation/comma_if.hpp>
38 #include <boost/preprocessor/arithmetic/add.hpp>
39 #include <boost/preprocessor/comparison/equal.hpp>
40 #include <boost/preprocessor/repetition/repeat.hpp>
41 #define DEF_PROPERTY_use_properties(clname, type, condition, name, getter) \
42  private:\
43  inline const type BOOST_PP_CAT(get_, name)() const { assert(condition); return getter }\
44  public:\
45  property<clname, const type> name;
46 #define DEF_PROPERTY_no_use_properties(clname, type, condition, name, getter) type name = NotAssigned;
47 #define DEF_PROPERTY_impl(clname, type, condition, name, getter) \
48  BOOST_PP_IF(BOOST_PP_CAT(clname, _USE_PROPERTIES), DEF_PROPERTY_use_properties, DEF_PROPERTY_no_use_properties)(clname, type, condition, name, getter)
49 #define DEF_PROPERTY(r,data,impl) \
50  BOOST_PP_IF(BOOST_PP_EQUAL(BOOST_PP_SEQ_SIZE(impl), 2), DEF_PROPERTY_impl, CONDITIONAL_PROPERTY_impl)(\
51  BOOST_PP_SEQ_ELEM(0, data), BOOST_PP_SEQ_ELEM(1, data), BOOST_PP_SEQ_ELEM(2, data),\
52  BOOST_PP_SEQ_ELEM(BOOST_PP_ADD(0, BOOST_PP_EQUAL(BOOST_PP_SEQ_SIZE(impl), 3)), impl),\
53  BOOST_PP_SEQ_ELEM(BOOST_PP_ADD(1, BOOST_PP_EQUAL(BOOST_PP_SEQ_SIZE(impl), 3)), impl)\
54  )
55 #define DEF_PROPERTY2(z,i,dataprops) DEF_PROPERTY2_impl(BOOST_PP_SEQ_ELEM(0,dataprops), BOOST_PP_SEQ_ELEM(i,BOOST_PP_SEQ_ELEM(1,dataprops)))
56 #define DEF_PROPERTY2_impl(data, impl) BOOST_PP_IF(BOOST_PP_CAT(BOOST_PP_SEQ_ELEM(0, data), _USE_PROPERTIES), DEF_PROPERTY_use_properties, DEF_PROPERTY_no_use_properties)(\
57  BOOST_PP_SEQ_ELEM(0, data), BOOST_PP_SEQ_ELEM(1, data), BOOST_PP_SEQ_ELEM(2, data),\
58  BOOST_PP_SEQ_ELEM(0, impl), BOOST_PP_SEQ_ELEM(1, impl)\
59  )
60 #define CONDITIONAL_PROPERTY_impl(clname, type, old_conds, condition, props) BOOST_PP_REPEAT(BOOST_PP_SEQ_SIZE(props), DEF_PROPERTY2, ((clname)(type)((old_conds) && (condition)))(props))
61 #define CONDITIONAL_PROPERTIES(condition, props) (1)(condition)(props)
62 
63 #define ASSERT_INITIALIZED(r,data,prop) assert(BOOST_PP_SEQ_ELEM(0,prop) != NotAssigned);
64 #define PROP_INIT_impl(check_these, pname, getter) BOOST_PP_SEQ_FOR_EACH(ASSERT_INITIALIZED, ~, check_these) pname = getter
65 #define PROP_INIT_cond_impl(check_these, cond, props) if(cond) { BOOST_PP_SEQ_FOR_EACH(PROP_INIT2, ~, props) }
66 
67 #define PROP_INIT2(r,data,impl) BOOST_PP_SEQ_ELEM(0, impl) = BOOST_PP_SEQ_ELEM(1, impl)
68 
69 #define PROP_INIT(r,data,i,impl) BOOST_PP_IF(BOOST_PP_EQUAL(BOOST_PP_SEQ_SIZE(impl), 2), PROP_INIT_impl, PROP_INIT_cond_impl)(\
70  BOOST_PP_SEQ_FIRST_N(i,data), \
71  BOOST_PP_SEQ_ELEM(BOOST_PP_ADD(0, BOOST_PP_EQUAL(BOOST_PP_SEQ_SIZE(impl), 3)), impl),\
72  BOOST_PP_SEQ_ELEM(BOOST_PP_ADD(1, BOOST_PP_EQUAL(BOOST_PP_SEQ_SIZE(impl), 3)), impl)\
73 )
74 
75 #define PROPS_ASSIGN_nocond(clname, prop) BOOST_PP_SEQ_ELEM(0,prop).set_target_getter_setter(this, BOOST_PP_CAT(&clname::get_, BOOST_PP_SEQ_ELEM(0,prop)));
76 #define PROPS_ASSIGN_cond_impl(r, i, prop) BOOST_PP_SEQ_ELEM(0,prop).set_target(this);
77 #define PROPS_ASSIGN_cond(clname, prop) BOOST_PP_SEQ_FOR_EACH(PROPS_ASSIGN_cond_impl, ~, BOOST_PP_SEQ_ELEM(2,prop))
78 #define PROPS_ASSIGN_impl(r,data,i,prop) BOOST_PP_IF(BOOST_PP_EQUAL(BOOST_PP_SEQ_SIZE(prop), 2), PROPS_ASSIGN_nocond, PROPS_ASSIGN_cond)(data, prop)
79 #define PROPS_ASSIGN(clname, props) BOOST_PP_SEQ_FOR_EACH_I(PROPS_ASSIGN_impl, clname, props)
80 #define NO_PROPS_ASSIGN(clname, props) init();
81 #define DEFINE_PROPERTIES(clname, type, props, init_function_name) \
82  /* Create the member variables */ \
83  BOOST_PP_SEQ_FOR_EACH_R(1, DEF_PROPERTY, (clname)(type)(true), props) \
84  /* create the init() function */ \
85  BOOST_PP_IF(BOOST_PP_CAT(clname, _USE_PROPERTIES), \
86  private: \
87  void init_function_name() { \
88  BOOST_PP_SEQ_FOR_EACH_I(PROPS_ASSIGN_impl, clname, props) \
89  } \
90  public: \
91  , \
92  private: \
93  void init_function_name() { \
94  BOOST_PP_SEQ_FOR_EACH_I( PROP_INIT, props, props) \
95  } \
96  public: \
97  )
98 
99 #define ShellData_USE_PROPERTIES 1
100 struct ShellData : public BasisElementData {
101 
102 
103  DEFINE_PROPERTIES(ShellData, int,
104  ((bfoff)(
105  basis->shell_to_function(index);
106  ))
107  ((nbf)(
108  basis->shell(index).nfunction();
109  ))
110  ((center)(
111  basis->shell_to_center(index);
112  ))
113  ((atom_shoff)(
114  basis->shell_on_center(center, 0);
115  ))
116  ((atom_bfoff)(
117  basis->shell_to_function(atom_shoff);
118  ))
119  ((atom_nsh)(
120  basis->nshell_on_center(center);
121  ))
122  ((atom_nbf)(
123  basis->nbasis_on_center(center);
124  ))
125  ((shoff_in_atom)(
126  index - atom_shoff;
127  ))
128  ((bfoff_in_atom)(
129  bfoff - atom_bfoff;
130  ))
131  ((atom_last_function)(
132  atom_bfoff + atom_nbf - 1;
133  ))
134  ((atom_last_shell)(
135  atom_shoff + atom_nsh - 1;
136  ))
137  ((last_function)(
138  bfoff + nbf - 1;
139  ))
140  (CONDITIONAL_PROPERTIES(dfbasis != 0,
141  ((atom_dfshoff)(
142  dfbasis->shell_on_center(center, 0);
143  ))
144  ((atom_dfbfoff)(
145  dfbasis->shell_to_function(atom_dfshoff);
146  ))
147  ((atom_dfnsh)(
148  dfbasis->nshell_on_center(center);
149  ))
150  ((atom_dfnbf)(
151  dfbasis->nbasis_on_center(center);
152  ))
153  ((atom_df_last_function)(
154  atom_dfbfoff + atom_dfnbf - 1; ))
155  ((atom_df_last_shell)(
156  atom_dfshoff + atom_dfnsh - 1;
157  ))
158  )),
159  init
160  );
161 
162  ShellData()
163  : BasisElementData(NotAssigned, 0, 0)
164  { }
165 
166  ShellData(
167  int idx,
168  GaussianBasisSet* basis,
169  GaussianBasisSet* dfbasis = 0
170  ) : BasisElementData(idx, basis, dfbasis)
171  {
172  init();
173  }
174 
175  ShellData(
176  const ShellData& other
177  ) : BasisElementData(other.index, other.basis, other.dfbasis)
178  {
179  init();
180  }
181 
182  const ShellData& operator++()
183  {
184  ++index;
185  #if !ShellData_USE_PROPERTIES
186  init();
187  #endif
188  return *this;
189  }
190 
191  void set_index(int idx) {
192  index = idx;
193  #if !ShellData_USE_PROPERTIES
194  init();
195  #endif
196  }
197 
198  const ShellData& operator*() const { return *this; }
199 
200  ShellData&
201  operator=(const ShellData& other)
202  {
203  index = other.index;
204  basis = other.basis;
205  dfbasis = other.dfbasis;
206  if(&other != this){
207  init();
208  }
209  return *this;
210  }
211 
212  operator int() { ASSERT_SHELL_BOUNDS; return index; }
213  operator const int() const { ASSERT_SHELL_BOUNDS; return index; }
214 
215 };
216 
217 //============================================================================//
218 
219 namespace detail {
220 
221 template <typename DataContainer, typename Iterable>
223 
224  protected:
225 
226  enum { NoLastIndex = -1, NoFirstIndex = -2 };
227 
228  GaussianBasisSet* basis_;
229  GaussianBasisSet* dfbasis_;
230  const Iterable& indices_;
231 
232  private:
233 
234  std::shared_ptr<Iterable> created_index_iterator_;
235 
236  public:
237 
239  typedef decltype(Iterable().begin()) index_iterator_type;
241 
242  //basis_iterable() : basis_(0), dfbasis_(0), indices_() { }
243 
245  GaussianBasisSet* basis,
246  GaussianBasisSet* dfbasis, //OptionalRefParameter<GaussianBasisSet> dfbasis,
247  int first_index,
248  int last_index
249  ) : basis_(basis), dfbasis_(dfbasis),
250  created_index_iterator_(std::make_shared<Iterable>(
251  boost::counting_range(
252  first_index, last_index == NoLastIndex ? DataContainer::max_index(basis_) : last_index
253  )
254  )),
255  indices_(*created_index_iterator_)
256  { }
257 
259  const Ref<GaussianBasisSet>& basis,
260  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
261  ) : basis_iterable(basis, dfbasis, 0, NoLastIndex)
262  { }
263 
264  template<typename iterable_like, typename=typename std::enable_if<
265  std::is_base_of<Iterable, typename std::decay<iterable_like>::type>::value
266  >::type
267  >
269  const Ref<GaussianBasisSet>& basis,
270  const OptionalRefParameter<GaussianBasisSet>& dfbasis,
271  const iterable_like& indices
272  ) : basis_(basis), dfbasis_(dfbasis), indices_(indices)
273  { }
274 
275  template<typename iterable_like, typename=typename std::enable_if<
276  std::is_base_of<Iterable, typename std::decay<iterable_like>::type>::value
277  >::type
278  >
280  const Ref<GaussianBasisSet>& basis,
281  const iterable_like& indices,
282  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
283  ) : basis_(basis), dfbasis_(dfbasis), indices_(indices)
284  { }
285 
287  const Ref<GaussianBasisSet>& basis,
288  int last_index,
289  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
290  ) : basis_iterable(basis, dfbasis, 0, last_index)
291  { }
292 
293  template<typename int_like, typename=typename std::enable_if<
294  std::is_convertible<int_like, int>::value
295  >::type
296  >
298  const Ref<GaussianBasisSet>& basis,
299  const Ref<GaussianBasisSet>& dfbasis,
300  int_like last_index
301  ) : basis_iterable(basis, dfbasis, 0, (int)last_index)
302  { }
303 
305  const Ref<GaussianBasisSet>& basis,
306  int first_index,
307  int last_index,
308  const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
309  ) : basis_iterable(basis, dfbasis, first_index, last_index)
310  { }
311 
313  const DataContainer& first,
314  const DataContainer& last
315  ) : basis_iterable(first.basis, first.dfbasis, Iterable(first, last))
316  { }
317 
318 
319  value_type begin() const;
320  value_type last() const;
321  value_type end() const;
322 
323 };
324 
325 } // end namespace detail
326 
327 //============================================================================//
328 
329 template <typename Iterable>
330 class shell_iterable : public detail::basis_iterable<ShellData, Iterable> {
331 
332  public:
333 
336  typedef decltype(Iterable().begin()) index_iterator_type;
338 
339  // Not sure if I need these anymore or not...
340  shell_iterable(self_type&) = default;
341  shell_iterable(const self_type&) = default;
342 
343  explicit shell_iterable(const ShellBlockIterator<Iterable>& block);
344 
346 
347 };
348 
353 template<typename ShellIncrementable>
355 make_shell_iterable(ShellIncrementable first, ShellIncrementable last) {
356 
357  static_assert(not std::is_convertible<ShellIncrementable, BasisFunctionData>::value,
358  "Can't construct a shell_iterable using BasisFunctionData objects."
359  );
360 
362  first.basis,
363  first.dfbasis,
364  boost::iterator_range<boost::counting_iterator<ShellIncrementable>>(
365  first, last
366  )
367  );
368 }
369 
371 
372 //============================================================================//
373 
374 template <typename Iterable = int_range>
375 class function_iterable : public detail::basis_iterable<BasisFunctionData, Iterable> {
376 
377  static constexpr int NotAssigned = -1;
378 
379  int block_offset = NotAssigned;
380 
381  public:
382 
385  typedef decltype(Iterable().begin()) index_iterator_type;
387 
388  // Not sure if I need these anymore or not...
389  function_iterable() = default;
391  function_iterable(const function_iterable&) = default;
392 
393  explicit function_iterable(const ShellData&);
394 
396 
398 
399 };
400 
402 
403 //============================================================================//
404 
405 template<typename Iterable>
406 class shell_block_iterable : public detail::basis_iterable<ShellBlockIterator<Iterable>, Iterable> {
407 
408  int target_size = DEFAULT_TARGET_BLOCK_SIZE;
409 
410  // Composed using bitwise or of BlockCompositionRequirement enums
411  int reqs = SameCenter;
412 
413  public:
414 
417  typedef decltype(Iterable().begin()) index_iterator_type;
419 
420  shell_block_iterable() = default;
421  shell_block_iterable(self_type&) = default;
422  shell_block_iterable(const self_type&) = default;
423 
424  using detail::basis_iterable<ShellBlockIterator<Iterable>, Iterable>::basis_iterable;
425 
426  const shell_block_iterable& requiring(int in_reqs) {
427  reqs = in_reqs;
428  return *this;
429  }
430 
431  const shell_block_iterable& with_target_size(int new_target_size) {
432  target_size = new_target_size;
433  return *this;
434  }
435 
436  value_type begin() const;
437  value_type last() const;
438  value_type end() const;
439 
440 };
441 
443 
444 template<typename Iterator, typename ParentContainer>
445 class IterableBasisElementData : public ParentContainer {
446 
447  Iterator spot;
448 
449  public:
450 
452  typedef boost::forward_traversal_tag iterator_category;
453  typedef typename Iterator::difference_type difference_type;
454 
455  template<typename... Args>
457  Iterator index_iter,
458  GaussianBasisSet* basis,
459  GaussianBasisSet* dfbasis,
460  Args&&... args
461  ) : spot(index_iter), ParentContainer(*index_iter, basis, dfbasis, std::forward<Args>(args)...)
462  { }
463 
464  const self_type& operator++()
465  {
466  ++spot;
467  ParentContainer::index = *spot;
468  ParentContainer::init();
469  return *this;
470  }
471 
472  const self_type& operator*() const {
473  return *this;
474  }
475 
476 };
477 
478 template<typename range_type>
480 
481  public:
482 
483  static constexpr int NoMaximumBlockSize = -20;
484 
485  // Forward declaration of inner class
486  class Skeleton;
487 
488  // Construct a ShellBlockIterator that spans an entire basis
490  GaussianBasisSet* basis,
491  GaussianBasisSet* dfbasis = 0
492  ) : BasisElementData(NotAssigned, basis, dfbasis),
493  possible_shell_iter(basis, dfbasis, first_shell, basis->nshell()),
494  target_size(NoMaximumBlockSize),
495  reqs(NoRestrictions),
496  first_shell(possible_shell_iter.begin()),
497  last_shell(possible_shell_iter.end()), // for now; init will fix this
498  shell_iter(first_shell, last_shell) // for now; init() will fix this
499  {
500  init();
501  }
502 
504  int first_index,
505  GaussianBasisSet* basis,
506  GaussianBasisSet* dfbasis,
507  int reqs = SameCenter,
508  int target_size = DEFAULT_TARGET_BLOCK_SIZE
509  ) : BasisElementData(NotAssigned, basis, dfbasis),
510  possible_shell_iter(basis, dfbasis, first_index, basis->nshell()),
511  target_size(target_size),
512  reqs(reqs),
513  first_shell(possible_shell_iter.begin()),
514  last_shell(possible_shell_iter.end()), // for now; init will fix this
515  shell_iter(first_shell, last_shell) // for now; init() will fix this
516  {
517  init();
518  }
519 
520  template<typename iterable_like, typename=typename std::enable_if<
521  std::is_base_of<Iterable, typename std::decay<iterable_like>::type>::value
522  && std::is_convertible<iterable_like, Iterable>::value
523  >::type>
525  const iterable_like& iterable,
526  GaussianBasisSet* basis,
527  GaussianBasisSet* dfbasis,
528  int reqs = SameCenter,
529  int target_size = DEFAULT_TARGET_BLOCK_SIZE
530  ) : BasisElementData(NotAssigned, basis, dfbasis),
531  possible_shell_iter(basis, dfbasis, iterable),
532  target_size(target_size),
533  reqs(reqs),
534  first_shell(possible_shell_iter.begin()),
535  last_shell(possible_shell_iter.end()), // for now; init() will fix this
536  shell_iter(first_shell, last_shell) // for now; init() will fix this
537  {
538  init();
539  }
540 
541 
543 
550  {
551  return first_shell.index != other.first_shell.index or nshell != other.nshell;
552  }
553 
554  const ShellBlockIterator& operator*() const { return *this; }
555 
556  GaussianBasisSet* basis;
557  GaussianBasisSet* dfbasis;
558  shell_iterable<Iterable> shell_iter;
559  typename shell_iterable<Iterable>::value_type first_shell;
560  typename shell_iterable<Iterable>::value_type last_shell;
561 
562  int nbf;
563  int bfoff;
564  int nshell;
565  int last_function;
566  int center = -1;
567  int atom_bfoff = -1;
568  int atom_shoff = -1;
569  int atom_nsh = -1;
570  int atom_nbf = -1;
571  int bfoff_in_atom = -1;
572  int shoff_in_atom = -1;
573  int atom_last_function = -1;
574  int atom_last_shell = -1;
575  int atom_dfshoff = -1;
576  int atom_dfbfoff = -1;
577  int atom_dfnbf = -1;
578  int atom_dfnsh = -1;
579  int atom_df_last_function = -1;
580  int atom_df_last_shell = -1;
581 
582  static int max_index(GaussianBasisSet* basis) { return basis->nshell(); }
583 
584  protected:
585 
586  shell_iterable<Iterable> possible_shell_iter;
587  int first_index;
588  int target_size;
589  int reqs;
590 
591  void init();
592 
593  public:
594 
595  class Skeleton {
596 
597  public:
598 
599  int first_index;
600  int size;
601  int reqs;
602  GaussianBasisSet* basis;
603  GaussianBasisSet* dfbasis;
604 
605  Skeleton(
606  const ShellBlockIterator& data
607  ) : first_index(data.first_index),
608  size(data.nbf),
609  basis(data.basis),
610  dfbasis(data.dfbasis),
611  reqs(data.reqs)
612  { }
613 
614  bool operator<(const ShellBlockIterator::Skeleton& other) const {
615  return first_index < other.first_index
616  or (first_index == other.first_index and size < other.size);
617  }
618 
619  bool operator!=(const ShellBlockIterator::Skeleton& other) const {
620  return first_index != other.first_index or size != other.size;
621  }
622 
623  };
624 
625 };
626 
628 
629 template<typename Iterator, typename Iterable>
630 class IterableBasisElementData<Iterator, ShellBlockIterator<Iterable>> : public ShellBlockIterator<Iterable> {
631 
634 
646 
647  protected:
648 
649  Iterator spot;
650 
651  public:
652 
653  // TODO avoid nesting of boost::counting_range types?
654  template<typename... Args>
656  const Iterable& index_iter,
657  GaussianBasisSet* basis,
658  GaussianBasisSet* dfbasis,
659  Args&&... args
660  ) : spot(index_iter.begin()),
662  index_iter,
663  basis, dfbasis,
664  std::forward<Args>(args)...
665  )
666  { }
667 
668  const self_type& operator++();
669 
670  const ShellBlockIterator<Iterable>& operator*() const { return *this; }
671 
672  const self_type& advance_to_last_block();
673 
674 };
675 
676 template <typename DataContainer, typename Iterable>
679 {
680  return value_type(indices_.begin(), basis_, dfbasis_);
681 }
682 
683 template <typename DataContainer, typename Iterable>
686 {
687  auto last = indices_.end();
688  --last;
689  return value_type(last, basis_, dfbasis_);
690 }
691 
692 template <typename DataContainer, typename Iterable>
695 {
696  return value_type(indices_.end(), basis_, dfbasis_);
697 }
698 
699 template <typename Iterable>
702 {
703  return value_type(
704  base_type::indices_,
705  base_type::basis_, base_type::dfbasis_,
706  reqs, target_size
707  );
708 }
709 
710 template <typename Iterable>
713 {
714  return value_type(
715  Iterable(base_type::indices_.end(), base_type::indices_.end()),
716  base_type::basis_, base_type::dfbasis_,
717  reqs, target_size
718  );
719 }
720 
721 template <typename Iterable>
724 {
725  value_type rv(
726  base_type::indices_.begin(),
727  base_type::indices_.end(),
728  base_type::basis_, base_type::dfbasis_,
729  reqs, target_size
730  );
731  rv.advance_to_last_block();
732  return rv;
733 }
734 
735 template<typename Iterator, typename Iterable>
738 {
739 
740  //----------------------------------------//
741  // Treat the end iterator specially
742  if(possible_shell_iter.begin() == possible_shell_iter.end()){
743  throw ProgrammingError(
744  "Already at end; can't advance to last block, which is before end.",
745  __FILE__, __LINE__
746  );
747  }
748  auto iter_tmp = possible_shell_iter;
749  const auto& last_possible_shell = possible_shell_iter.last();
750  //----------------------------------------//
751  while(iter_tmp.begin() != iter_tmp.end()) {
752  first_shell = iter_tmp.begin();
753  nbf = 0;
754  nshell = 0;
755  //----------------------------------------//
756  const int first_center = first_shell.center;
757  auto& first_am = basis->shell(first_shell).am();
758  for(auto ish : possible_shell_iter){
759  if(
760  // Same center condition
761  ((reqs & SameCenter) and ish.center != first_center)
762  or
763  // Same angular momentum condition
764  ((reqs & SameAngularMomentum) and
765  basis->shell(ish).am() != first_am) or
766  // Maximum block size overflow condition
767  (target_size != NoMaximumBlockSize and nbf >= target_size)
768  ){
769  // Should always have at least one shell, unless I mess up
770  // the conditional expression above. Assert this fact:
771  assert(nshell);
772  // We're done making the block. `ish` belongs to the next block.
773  break;
774  }
775  else{
776  // This could be the last shell. Overwrite the last shell member
777  // everytime we get here
778  last_shell = ish;
779  ++nshell;
780  nbf += ish.nbf;
781  }
782  }
783  auto new_first_shell = last_shell;
784  new_first_shell++;
785  iter_tmp = make_shell_iterable(new_first_shell, last_possible_shell);
786  }
787  //----------------------------------------//
788  possible_shell_iter = make_shell_iterable(first_shell, last_possible_shell);
789  init();
790 }
791 
792 template<typename Iterable>
793 inline
795  const typename ShellBlockIterator<Iterable>::Skeleton& sk
796 ) : ShellBlockIterator<Iterable>(sk.first_index, sk.basis, sk.dfbasis, sk.reqs, sk.size)
797 { }
798 
799 
800 template<typename Iterator, typename Iterable>
803 {
804  spot += nshell;
805  possible_shell_iter = make_shell_iter(spot, possible_shell_iter.last());
806  init();
807  return *this;
808 }
809 
811 // compute K attic
812 
813 // three body parts
814 // Compute B
815  /* Fail-safe explicit loops version
816  //for(auto sigma : function_range(obs)) {
817  // B_mus[mu.bfoff_in_shell].row(sigma) += 2.0 *
818  // D.row(sigma).segment(jsh.bfoff, jsh.nbf) *
819  // g3.middleRows(mu.bfoff_in_shell*jsh.nbf, jsh.nbf);
820  // //for(auto X : function_range(Xblk)) {
821  // // B_mus[mu.bfoff_in_shell](sigma, X.bfoff_in_block) += 2.0 *
822  // // D.row(sigma).segment(jsh.bfoff, jsh.nbf) *
823  // // g3.col(X.bfoff_in_block).segment(mu.bfoff_in_shell*jsh.nbf, jsh.nbf);
824  // // //for(auto rho : function_range(jsh)) {
825  // // // B_mus[mu.bfoff_in_shell](sigma, X.bfoff_in_block) += 2.0 *
826  // // // g3(mu.bfoff_in_shell*jsh.nbf + rho.bfoff_in_shell, X.bfoff_in_block)
827  // // // * D(rho, sigma);
828  // // //}
829  // //
830  // //}
831  //}
832  */
833 // K Contributeions
834 
835  /* Failsafe version:
836  //for(auto nu : iter_functions_on_center(obs, Y.center)) {
837  // //for(auto sigma : function_range(obs)) {
838  // // Kt_part(mu, nu) -= 0.5 * C_Y(nu.bfoff_in_atom, sigma) * B_mus[mu.bfoff_in_shell](sigma, Y.bfoff_in_shell);
839  // //}
840  // Kt_part(mu, nu) -= 0.5 * C_Y.row(nu.bfoff_in_atom) * B_mus[mu.bfoff_in_shell].col(Y.bfoff_in_shell);
841  //}
842  //----------------------------------------//
843  //for(auto nu : function_range(obs)) {
844  // if(nu.center != Y.center) {
845  // //for(auto sigma : iter_functions_on_center(obs, Y.center)) {
846  // // Kt_part(mu, nu) -= 0.5 * C_Y(sigma.bfoff_in_atom, nu) * B_mus[mu.bfoff_in_shell](sigma, Y.bfoff_in_shell);
847  // //}
848  // Kt_part(mu, nu) -= 0.5 * C_Y.col(nu).transpose() * B_mus[mu.bfoff_in_shell].col(Y.bfoff_in_shell).segment(
849  // obs->shell_to_function(obs->shell_on_center(Y.center, 0)),
850  // obs->nbasis_on_center(Y.center)
851  // );
852  // }
853  //}
854  */
855  //----------------------------------------//
856  #if OLD_K
857  for(auto ksh : shell_range(obs, dfbs_)) {
858  for(auto lsh : iter_significant_partners(ksh)) {
859  if(ksh.center != Xblk.center && lsh.center != Xblk.center) continue;
860  for(auto nu : function_range(ksh)) {
861  for(auto sigma : function_range(lsh)) {
862  CoefContainer c_ptr;
863  if(lsh <= ksh) {
864  IntPair nu_sigma(nu, sigma);
865  assert(coefs_.find(nu_sigma) != coefs_.end());
866  auto coef_pair = coefs_[nu_sigma];
867  c_ptr = ksh.center == Xblk.center ? coef_pair.first : coef_pair.second;
868  }
869  else {
870  IntPair nu_sigma(sigma, nu);
871  assert(coefs_.find(nu_sigma) != coefs_.end());
872  auto coef_pair = coefs_[nu_sigma];
873  c_ptr = lsh.center == Xblk.center ? coef_pair.first : coef_pair.second;
874  }
875  auto& C = *c_ptr;
876  for(auto X : function_range(Xblk)) {
877  for(auto mu : function_range(ish)) {
878  Kt_part(mu, nu) += C[X.bfoff_in_atom] * B_mus[mu.bfoff_in_shell](sigma, X.bfoff_in_block);
879  } // end loop over mu in ish
880  } // end loop over X in Xblk
881  } // end loop over sigma in lsh
882  } // end loop over nu in ksh
883  } // end loop over lsh
884  } // end loop over ksh
885  #endif
886 
887 // Full OLD_K section
888  #if OLD_K
889  ShellData jsh;
890  while(get_shell_pair(ish, jsh, SignificantPairs)){
891  const double pf = 2.0;
892  /*-----------------------------------------------------*/
893  /* Setup {{{2 */ #if 2 // begin fold
894  std::vector<Eigen::MatrixXd> A_mus(ish.nbf);
895  std::vector<Eigen::MatrixXd> A_rhos((ish == jsh) ? 0 : (int)jsh.nbf);
896  std::vector<Eigen::MatrixXd> B_mus(ish.nbf);
897  std::vector<Eigen::MatrixXd> B_rhos((ish == jsh) ? 0 : (int)jsh.nbf);
898  std::vector<Eigen::MatrixXd> dt_mus(ish.nbf);
899  std::vector<Eigen::MatrixXd> dt_rhos((ish == jsh) ? 0 : (int)jsh.nbf);
900  std::vector<Eigen::MatrixXd> gt_mus(ish.nbf);
901  std::vector<Eigen::MatrixXd> gt_rhos((ish == jsh) ? 0 : (int)jsh.nbf);
902  for(auto mu : function_range(ish)){
903  A_mus[mu.bfoff_in_shell].resize(nbf, dfnbf);
904  A_mus[mu.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
905  B_mus[mu.bfoff_in_shell].resize(nbf, dfnbf);
906  B_mus[mu.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
907  dt_mus[mu.bfoff_in_shell].resize(nbf, dfnbf);
908  dt_mus[mu.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
909  gt_mus[mu.bfoff_in_shell].resize(nbf, dfnbf);
910  gt_mus[mu.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
911  }
912  if(ish != jsh){
913  for(auto rho : function_range(jsh)){
914  A_rhos[rho.bfoff_in_shell].resize(nbf, dfnbf);
915  A_rhos[rho.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
916  B_rhos[rho.bfoff_in_shell].resize(nbf, dfnbf);
917  B_rhos[rho.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
918  dt_rhos[rho.bfoff_in_shell].resize(nbf, dfnbf);
919  dt_rhos[rho.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
920  gt_rhos[rho.bfoff_in_shell].resize(nbf, dfnbf);
921  gt_rhos[rho.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, dfnbf);
922  }
923  }
924  /*******************************************************/ #endif //2}}}
925  /*-----------------------------------------------------*/
926  /* Compute B for functions in ish, jsh {{{2 */ #if 2 // begin fold
927  if(ithr == 0) timer.enter("build B");
928  for(int kshdf = 0; kshdf < dfbs_->nshell(); ++kshdf){
929  const int kshdf_bfoff = dfbs_->shell_to_function(kshdf);
930  const int kshdf_nbf = dfbs_->shell(kshdf).nfunction();
931  //----------------------------------------//
932  // Get the integrals for ish, jsh, ksh
933  if(ithr == 0) timer.enter("compute ints");
934  std::shared_ptr<Eigen::MatrixXd> g_part = ints_to_eigen(
935  ish, jsh, kshdf,
936  eris_3c_[ithr],
937  coulomb_oper_type_
938  );
939  //----------------------------------------//
940  // B_mus
941  if(ithr == 0) timer.change("B_mus");
942  for(int ibf = 0; ibf < ish.nbf; ++ibf){
943  B_mus[ibf].middleCols(kshdf_bfoff, kshdf_nbf) +=
944  pf * D.middleCols(jsh.bfoff, jsh.nbf)
945  * g_part->middleRows(ibf * jsh.nbf, jsh.nbf);
946  } // end loop over mu
947  //----------------------------------------//
948  if(ithr == 0) timer.change("B_rhos");
949  if(ish != jsh) {
950  // Do the contribution to the other way around
951  for(int ibf = 0, mu = ish.bfoff; ibf < ish.nbf; ++ibf, ++mu){
952  for(int jbf = 0; jbf < jsh.nbf; ++jbf){
953  //----------------------------------------//
954  B_rhos[jbf].middleCols(kshdf_bfoff, kshdf_nbf) +=
955  pf * D.col(mu) * g_part->row(ibf * jsh.nbf + jbf);
956  //----------------------------------------//
957  }
958  }
959  } // end if ish != jsh
960  //----------------------------------------//
961  if(ithr == 0) timer.exit();
962  } // end loop over kshdf
963  //----------------------------------------//
964  if(xml_debug) {
965  for(auto mu : function_range(ish)){
966  write_as_xml("B_part", B_mus[mu.bfoff_in_shell], std::map<std::string, int>{ {"mu", mu}, {"jsh", jsh} } );
967  }
968  if(ish != jsh){
969  for(auto rho : function_range(jsh)){
970  write_as_xml("B_part", B_rhos[rho.bfoff_in_shell], std::map<std::string, int>{ {"mu", rho}, {"jsh", jsh} } );
971  }
972  }
973  }
974  //----------------------------------------//
975  if(ithr == 0) timer.exit("build B");
976  /*******************************************************/ #endif //2}}}
977  /*-----------------------------------------------------*/
978  /* Compute d_tilde and g_tilde for bfs in ish,jsh {{{2 */ #if 2 // begin fold
979  //if(ithr == 0) timer.enter("build d_tilde");
980  //for(auto mu : function_range(ish)){
981  // //----------------------------------------//
982  // // Compute d_tilde for a given mu
983  // for(auto rho : function_range(jsh)){
984  // //----------------------------------------//
985  // // Get the coefficients
986  // IntPair mu_rho(mu, rho);
987  // auto cpair = coefs_[mu_rho];
988  // VectorMap& Ca = *(cpair.first);
989  // VectorMap& Cb = *(cpair.second);
990  // //----------------------------------------//
991  // // dt_mus
992  // // column of D times Ca as row vector
993  // dt_mus[mu.bfoff_in_shell].middleCols(ish.atom_dfbfoff, ish.atom_dfnbf) +=
994  // pf * D.row(rho).transpose() * Ca.transpose();
995  // if(ish.center != jsh.center) {
996  // dt_mus[mu.bfoff_in_shell].middleCols(jsh.atom_dfbfoff, jsh.atom_dfnbf) +=
997  // pf * D.row(rho).transpose() * Cb.transpose();
998  // }
999  // //----------------------------------------//
1000  // if(ish != jsh) {
1001  // // dt_rhos
1002  // // same thing, utilizing C_{mu,rho} = C_{rho,mu} (accounting for ordering, etc.)
1003  // dt_rhos[rho.bfoff_in_shell].middleCols(ish.atom_dfbfoff, ish.atom_dfnbf) +=
1004  // pf * D.row(mu).transpose() * Ca.transpose();
1005  // if(ish.center != jsh.center) {
1006  // dt_rhos[rho.bfoff_in_shell].middleCols(jsh.atom_dfbfoff, jsh.atom_dfnbf) +=
1007  // pf * D.row(mu).transpose() * Cb.transpose();
1008  // }
1009  // }
1010  // //----------------------------------------//
1011  // } // end loop over rho
1012  // //----------------------------------------//
1013  //} // end loop over mu
1014  //if(ithr == 0) timer.exit("build d_tilde");
1015  /*******************************************************/ #endif //2}}}
1016  /*-----------------------------------------------------*/
1017  /* Get g_tilde contribution for bfs in ish,jsh {{{2 */ #if 2 // begin fold
1018  //if(ithr == 0) timer.enter("build g");
1019  /* TODO optimally, this needs to compute blocks of several shells at once and contract these blocks
1020  for(auto Ysh : shell_range(dfbs_)) {
1021  for(auto Xsh : shell_range(dfbs_, 0, Ysh)) {
1022  if(ithr == 0) timer.enter("compute ints");
1023  const double dfpf = Xsh == Ysh ? 1.0 : 2.0;
1024  Eigen::MatrixXd& g2_part = *ints_to_eigen(
1025  Xsh, Ysh,
1026  eris_2c_[ithr],
1027  coulomb_oper_type_
1028  );
1029  //----------------------------------------//
1030  if(ithr == 0) timer.change("gt_mu");
1031  for(auto mu : function_range(ish)) {
1032  gt_mus[mu.bfoff_in_shell].middleCols(Ysh.bfoff, Ysh.nbf).noalias() +=
1033  dfpf * dt_mus[mu.bfoff_in_shell].middleCols(Xsh.bfoff, Xsh.nbf) * g2_part;
1034  }
1035  //----------------------------------------//
1036  if(ish != jsh){
1037  if(ithr == 0) timer.change("gt_rho");
1038  for(auto rho : function_range(jsh)) {
1039  gt_rhos[rho.bfoff_in_shell].middleCols(Ysh.bfoff, Ysh.nbf).noalias() +=
1040  dfpf * dt_rhos[rho.bfoff_in_shell].middleCols(Xsh.bfoff, Xsh.nbf) * g2_part;
1041  } // end loop over rho
1042  }
1043  //----------------------------------------//
1044  if(ithr == 0) timer.exit();
1045  } // end loop over Xsh
1046  } // end loop over Ysh
1047  */
1048  //if(ithr == 0) timer.enter("gt_mu");
1049  //for(int mu = 0; mu < ish.nbf; ++mu){
1050  // //gt_mus[mu].noalias() += dt_mus[mu] * g2;
1051  // gt_mus[mu] = dt_mus[mu] * g2;
1052  //} // end loop over mu
1054  //if(ish != jsh){
1055  // if(ithr == 0) timer.change("gt_rho");
1056  // for(int rho = 0; rho < jsh.nbf; ++rho) {
1057  // //gt_rhos[rho].noalias() += dt_rhos[rho] * g2;
1058  // gt_rhos[rho] = dt_rhos[rho] * g2;
1059  // } // end loop over rho
1060  //}
1061  //----------------------------------------//
1062  // Subtract off the term three contribution to B_mus and B_rhos (result is called A in notes)
1063  //if(ithr == 0) timer.enter("compute A");
1064  //for(auto mu : function_range(ish)) {
1065  // A_mus[mu.bfoff_in_shell] = B_mus[mu.bfoff_in_shell] - 0.5 * dt_mus[mu.bfoff_in_shell] * g2; //gt_mus[mu.bfoff_in_shell];
1066  //}
1067  //if(ish != jsh) {
1068  // for(auto rho : function_range(jsh)) {
1069  // A_rhos[rho.bfoff_in_shell] = B_rhos[rho.bfoff_in_shell] - 0.5 * dt_rhos[rho.bfoff_in_shell] * g2; // gt_rhos[rho.bfoff_in_shell];
1070  // }
1071  //}
1072  //if(ithr == 0) timer.exit("compute A");
1073  //if(ithr == 0) timer.exit("build g");
1074  /*******************************************************/ #endif //2}}}
1075  /*-----------------------------------------------------*/
1076  /* Compute K contributions {{{2 */ #if 2 // begin fold
1077  if(ithr == 0) timer.enter("K contributions");
1078  //if(ithr == 0) timer.enter("misc");
1079  for(auto Y : function_range(dfbs_)) {
1080  Eigen::MatrixXd& C_Y = coefs_transpose_[Y];
1081  const int obs_atom_bfoff = obs->shell_to_function(obs->shell_on_center(Y.center, 0));
1082  const int obs_atom_nbf = obs->nbasis_on_center(Y.center);
1083  for(auto mu : function_range(ish)) {
1084  // B_mus[mu.bfoff_in_shell] is (nbf x Ysh.nbf)
1085  // C_Y is (Y.{obs_}atom_nbf x nbf)
1086  // result should be (Y.{obs_}atom_nbf x 1)
1087  Kt_part.row(mu).segment(obs_atom_bfoff, obs_atom_nbf).transpose() +=
1088  C_Y * B_mus[mu.bfoff_in_shell].col(Y);
1089  // The sigma <-> nu term
1090  Kt_part.row(mu).transpose() += C_Y.transpose()
1091  * B_mus[mu.bfoff_in_shell].col(Y).segment(obs_atom_bfoff, obs_atom_nbf);
1092  // Add back in the nu.center == Y.center part
1093  Kt_part.row(mu).segment(obs_atom_bfoff, obs_atom_nbf).transpose() -=
1094  C_Y.middleCols(obs_atom_bfoff, obs_atom_nbf).transpose()
1095  * B_mus[mu.bfoff_in_shell].col(Y).segment(obs_atom_bfoff, obs_atom_nbf);
1096  //----------------------------------------//
1097  /* Failsafe version:
1098  //for(auto nu : iter_functions_on_center(obs, Y.center)) {
1099  // //for(auto sigma : function_range(obs)) {
1100  // // Kt_part(mu, nu) -= 0.5 * C_Y(nu.bfoff_in_atom, sigma) * B_mus[mu.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1101  // //}
1102  // Kt_part(mu, nu) -= 0.5 * C_Y.row(nu.bfoff_in_atom) * B_mus[mu.bfoff_in_shell].col(Y.bfoff_in_shell);
1103  //}
1104  //----------------------------------------//
1105  //for(auto nu : function_range(obs)) {
1106  // if(nu.center != Y.center) {
1107  // //for(auto sigma : iter_functions_on_center(obs, Y.center)) {
1108  // // Kt_part(mu, nu) -= 0.5 * C_Y(sigma.bfoff_in_atom, nu) * B_mus[mu.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1109  // //}
1110  // Kt_part(mu, nu) -= 0.5 * C_Y.col(nu).transpose() * B_mus[mu.bfoff_in_shell].col(Y.bfoff_in_shell).segment(
1111  // obs->shell_to_function(obs->shell_on_center(Y.center, 0)),
1112  // obs->nbasis_on_center(Y.center)
1113  // );
1114  // }
1115  //}
1116  */
1117  //----------------------------------------//
1118  }
1119  if(ish != jsh){
1120  for(auto rho : function_range(jsh)) {
1121  // B_rhos[rho.bfoff_in_shell] is (nbf x Ysh.nbf)
1122  // C_Y is (Y.{obs_}atom_nbf x nbf)
1123  // result should be (Y.{obs_}atom_nbf x 1)
1124  Kt_part.row(rho).segment(obs_atom_bfoff, obs_atom_nbf).transpose() +=
1125  C_Y * B_rhos[rho.bfoff_in_shell].col(Y);
1126  // The sigma <-> nu term
1127  Kt_part.row(rho).transpose() += C_Y.transpose()
1128  * B_rhos[rho.bfoff_in_shell].col(Y).segment(obs_atom_bfoff, obs_atom_nbf);
1129  // Add back in the nu.center == Y.center part
1130  Kt_part.row(rho).segment(obs_atom_bfoff, obs_atom_nbf).transpose() -=
1131  C_Y.middleCols(obs_atom_bfoff, obs_atom_nbf).transpose()
1132  * B_rhos[rho.bfoff_in_shell].col(Y).segment(obs_atom_bfoff, obs_atom_nbf);
1133  //----------------------------------------//
1134  /* Failsafe version:
1135  //for(auto nu : function_range(obs)) {
1136  // for(auto sigma : function_range(obs)) {
1137  // if(nu.center == Y.center){
1138  // Kt_part(rho, nu) -= 0.5 * C_Y(nu.bfoff_in_atom, sigma) * dt_rhos[rho.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1139  // }
1140  // else if(sigma.center == Y.center){
1141  // Kt_part(rho, nu) -= 0.5 * C_Y(sigma.bfoff_in_atom, nu) * dt_rhos[rho.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1142  // }
1143  // }
1144  //}
1145  */
1146  //----------------------------------------//
1147  }
1148  }
1149  }
1150  /* Old version
1151  for(auto nu : function_range(obs, dfbs_)) {
1152  for(auto sigma : function_range(obsptr, dfbsptr, 0, nu)) {
1153  //----------------------------------------//
1154  // Get the coefficients
1155  IntPair nu_sigma(nu, sigma);
1156  auto cpair = coefs_[nu_sigma];
1157  VectorMap& Ca = *(cpair.first);
1158  VectorMap& Cb = *(cpair.second);
1159  //----------------------------------------//
1160  // Add in the contribution to Kt(mu, nu)
1161  for(auto mu : function_range(ish)) {
1162  Kt_part(mu, nu) +=
1163  B_mus[mu.bfoff_in_shell].row(sigma).segment(nu.atom_dfbfoff, nu.atom_dfnbf) * Ca;
1164  if(nu != sigma){
1165  Kt_part(mu, sigma) +=
1166  B_mus[mu.bfoff_in_shell].row(nu).segment(nu.atom_dfbfoff, nu.atom_dfnbf) * Ca;
1167  }
1168  if(nu.center != sigma.center) {
1169  Kt_part(mu, nu) +=
1170  B_mus[mu.bfoff_in_shell].row(sigma).segment(sigma.atom_dfbfoff, sigma.atom_dfnbf) * Cb;
1171  if(nu != sigma){
1172  Kt_part(mu, sigma) +=
1173  B_mus[mu.bfoff_in_shell].row(nu).segment(sigma.atom_dfbfoff, sigma.atom_dfnbf) * Cb;
1174  }
1175  }
1176  } // end loop over mu
1177  //----------------------------------------//
1178  // Add in the contribution to Kt(rho, nu)
1179  if(ish != jsh){
1180  for(auto rho : function_range(jsh)) {
1181  Kt_part(rho, nu) +=
1182  B_rhos[rho.bfoff_in_shell].row(sigma).segment(nu.atom_dfbfoff, nu.atom_dfnbf) * Ca;
1183  if(nu != sigma){
1184  Kt_part(rho, sigma) +=
1185  B_rhos[rho.bfoff_in_shell].row(nu).segment(nu.atom_dfbfoff, nu.atom_dfnbf) * Ca;
1186  }
1187  if(nu.center != sigma.center) {
1188  Kt_part(rho, nu) +=
1189  B_rhos[rho.bfoff_in_shell].row(sigma).segment(sigma.atom_dfbfoff, sigma.atom_dfnbf) * Cb;
1190  if(nu != sigma){
1191  Kt_part(rho, sigma) +=
1192  B_rhos[rho.bfoff_in_shell].row(nu).segment(sigma.atom_dfbfoff, sigma.atom_dfnbf) * Cb;
1193  }
1194  } // end if nu.center != sigma.center
1195  } // end loop over rho
1196  } // end if ish != jsh
1197  //----------------------------------------//
1198  } // end loop over sigma
1199  } // end loop over nu
1200  */
1201  if(ithr == 0) timer.exit("K contributions");
1202  /*******************************************************/ #endif //2}}}
1203  /*-----------------------------------------------------*/
1204  } // end while get_shelsh)
1205  #endif
1206 
1207 
1209 // Two body part of K
1210 
1211  //for(auto Y : function_range(Yblk)) {
1212  // Ct_mus[mu.bfoff_in_shell](rho, Y.bfoff_in_block) += 2.0 *
1213  // g2.row(Y.bfoff_in_block) * Cmu.segment(Xblk.bfoff_in_atom, Xblk.nbf);
1214  //}
1215  //for(auto X : function_range(Xblk)) {
1216  // Ct_mus[mu.bfoff_in_shell](rho, Y.bfoff_in_block) += 2.0 *
1217  // g2(Y.bfoff_in_block, X.bfoff_in_block) * Cmu(X.bfoff_in_atom);
1218  //}
1219 
1220  //for(auto Y : function_range(Yblk)) {
1221  // Ct_mus[mu.bfoff_in_shell](rho, Y.bfoff_in_block) += 2.0 *
1222  // g2.row(Y.bfoff_in_block) * Crho.segment(Xblk.bfoff_in_atom, Xblk.nbf);
1223  //}
1224  //for(auto X : function_range(Xblk)) {
1225  // Ct_mus[mu.bfoff_in_shell](rho, Y.bfoff_in_block) += 2.0 *
1226  // g2(Y.bfoff_in_block, X.bfoff_in_block) * Crho(X.bfoff_in_atom);
1227  //}
1228 
1229  #if OLD_K
1230  ShellData ish, jsh;
1231  while(get_shell_pair(ish, jsh, SignificantPairs)){
1232  for(auto Yblk : shell_block_range(dfbs_, 0, 0, NoLastIndex, NoRestrictions)) {
1233  //----------------------------------------//
1234  std::vector<Eigen::MatrixXd> Ct_mus(ish.nbf);
1235  for(auto mu : function_range(ish)) {
1236  Ct_mus[mu.bfoff_in_shell].resize(jsh.nbf, Yblk.nbf);
1237  Ct_mus[mu.bfoff_in_shell] = Eigen::MatrixXd::Zero(jsh.nbf, Yblk.nbf);
1238  }
1239  /*-----------------------------------------------------*/
1240  /* Form Ct_mu for each mu in ish {{{2 */ #if 2 // begin fold
1241  if(ithr == 0) timer.enter("form Ct_mu");
1242  //for(auto Xsh : iter_shells_on_center(dfbs_, ish.center)) {
1243  for(auto Xblk : iter_shell_blocks_on_center(dfbs_, ish.center)) {
1244  if(ithr == 0) timer.enter("compute ints");
1245  auto g2_part_ptr = ints_to_eigen(
1246  Yblk, Xblk,
1247  eris_2c_[ithr],
1248  coulomb_oper_type_
1249  );
1250  if(ithr == 0) timer.exit("compute ints");
1251  const auto& g2_part = *g2_part_ptr;
1252  for(auto mu : function_range(ish)){
1253  for(auto rho : function_range(jsh)){
1254  //----------------------------------------//
1255  // Get the coefficients
1256  IntPair mu_rho(mu, rho);
1257  auto cpair = coefs_[mu_rho];
1258  VectorMap& Ca = *(cpair.first);
1259  //----------------------------------------//
1260  Ct_mus[mu.bfoff_in_shell].row(rho.bfoff_in_shell).transpose() +=
1261  g2_part * Ca.segment(Xblk.bfoff_in_atom, Xblk.nbf);
1262  //----------------------------------------//
1263  } // end loop over functions rho in jsh
1264  } // end loop over functions mu in ish
1265  } // end loop of Xsh over shells on ish.center
1266  //----------------------------------------//
1267  if(ish.center != jsh.center){
1268  for(auto Xblk : iter_shell_blocks_on_center(dfbs_, jsh.center)) {
1269  auto g2_part_ptr = ints_to_eigen(
1270  Yblk, Xblk,
1271  eris_2c_[ithr],
1272  coulomb_oper_type_
1273  );
1274  const auto& g2_part = *g2_part_ptr;
1275  for(auto mu : function_range(ish)){
1276  for(auto rho : function_range(jsh)){
1277  //----------------------------------------//
1278  // Get the coefficients
1279  IntPair mu_rho(mu, rho);
1280  auto cpair = coefs_[mu_rho];
1281  VectorMap& Cb = *(cpair.second);
1282  //----------------------------------------//
1283  Ct_mus[mu.bfoff_in_shell].row(rho.bfoff_in_shell).transpose() +=
1284  g2_part * Cb.segment(Xblk.bfoff_in_atom, Xblk.nbf);
1285  //----------------------------------------//
1286  } // end loop over functions rho in jsh
1287  } // end loop over functions mu in ish
1288  } // end loop of Xsh over shells on jsh.center
1289  } // end if ish.center != jsh.center
1290  //----------------------------------------//
1291  /********************************************************/ #endif //2}}}
1292  /*-----------------------------------------------------*/
1293  /* Form dt_mu, dt_rho for functions in ish, jsh {{{2 */ #if 2 // begin fold
1294  if(ithr == 0) timer.change("form dt_mu, dt_rho");
1295  std::vector<Eigen::MatrixXd> dt_mus(ish.nbf);
1296  std::vector<Eigen::MatrixXd> dt_rhos((ish == jsh) ? 0 : (int)jsh.nbf);
1297  for(auto mu : function_range(ish)) {
1298  dt_mus[mu.bfoff_in_shell].resize(nbf, Yblk.nbf);
1299  dt_mus[mu.bfoff_in_shell] = 2.0 * D.middleCols(jsh.bfoff, jsh.nbf) * Ct_mus[mu.bfoff_in_shell];
1300  /*
1301  if(xml_debug) {
1302  Eigen::MatrixXd tmp(nbf, dfnbf);
1303  tmp = Eigen::MatrixXd::Zero(nbf, dfnbf);
1304  tmp.middleCols(Ysh.bfoff, Ysh.nbf) = dt_mus[mu.bfoff_in_shell];
1305  write_as_xml(
1306  "new_dt_part", tmp,
1307  std::map<std::string, int>{ {"mu", mu}, {"jsh", jsh} }
1308  );
1309  }
1310  */
1311  }
1312  if(ish != jsh){
1313  for(auto rho : function_range(jsh)) {
1314  dt_rhos[rho.bfoff_in_shell].resize(nbf, Yblk.nbf);
1315  dt_rhos[rho.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, Yblk.nbf);
1316  for(auto mu : function_range(ish)) {
1317  dt_rhos[rho.bfoff_in_shell] += 2.0 * D.col(mu) * Ct_mus[mu.bfoff_in_shell].row(rho.bfoff_in_shell);
1318  }
1319  /*
1320  if(xml_debug) {
1321  Eigen::MatrixXd tmp(nbf, dfnbf);
1322  tmp = Eigen::MatrixXd::Zero(nbf, dfnbf);
1323  tmp.middleCols(Ysh.bfoff, Ysh.nbf) = dt_rhos[rho.bfoff_in_shell];
1324  write_as_xml(
1325  "new_dt_part", tmp,
1326  std::map<std::string, int>{ {"mu", rho}, {"jsh", ish} }
1327  );
1328  }
1329  */
1330  }
1331  }
1332  /********************************************************/ #endif //2}}}
1333  /*-----------------------------------------------------*/
1334  /* Add contributions to Kt_part {{{2 */ #if 2 // begin fold
1335  if(ithr == 0) timer.change("K contributions");
1336  //----------------------------------------//
1337  for(auto Y : function_range(Yblk)) {
1338  Eigen::MatrixXd& C_Y = coefs_transpose_[Y];
1339  const int obs_atom_bfoff = obs->shell_to_function(obs->shell_on_center(Y.center, 0));
1340  const int obs_atom_nbf = obs->nbasis_on_center(Y.center);
1341  for(auto mu : function_range(ish)) {
1342  // dt_mus[mu.bfoff_in_shell] is (nbf x Ysh.nbf)
1343  // C_Y is (Y.{obs_}atom_nbf x nbf)
1344  // result should be (Y.{obs_}atom_nbf x 1)
1345  Kt_part.row(mu).segment(obs_atom_bfoff, obs_atom_nbf).transpose() -=
1346  0.5 * C_Y * dt_mus[mu.bfoff_in_shell].col(Y.bfoff_in_block);
1347  // The sigma <-> nu term
1348  Kt_part.row(mu).transpose() -= 0.5
1349  * C_Y.transpose()
1350  * dt_mus[mu.bfoff_in_shell].col(Y.bfoff_in_block).segment(obs_atom_bfoff, obs_atom_nbf);
1351  // Add back in the nu.center == Y.center part
1352  Kt_part.row(mu).segment(obs_atom_bfoff, obs_atom_nbf).transpose() += 0.5
1353  * C_Y.middleCols(obs_atom_bfoff, obs_atom_nbf).transpose()
1354  * dt_mus[mu.bfoff_in_shell].col(Y.bfoff_in_block).segment(obs_atom_bfoff, obs_atom_nbf);
1355  //----------------------------------------//
1356  /* Failsafe version:
1357  //for(auto nu : iter_functions_on_center(obs, Y.center)) {
1358  // //for(auto sigma : function_range(obs)) {
1359  // // Kt_part(mu, nu) -= 0.5 * C_Y(nu.bfoff_in_atom, sigma) * dt_mus[mu.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1360  // //}
1361  // Kt_part(mu, nu) -= 0.5 * C_Y.row(nu.bfoff_in_atom) * dt_mus[mu.bfoff_in_shell].col(Y.bfoff_in_shell);
1362  //}
1363  //----------------------------------------//
1364  //for(auto nu : function_range(obs)) {
1365  // if(nu.center != Y.center) {
1366  // //for(auto sigma : iter_functions_on_center(obs, Y.center)) {
1367  // // Kt_part(mu, nu) -= 0.5 * C_Y(sigma.bfoff_in_atom, nu) * dt_mus[mu.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1368  // //}
1369  // Kt_part(mu, nu) -= 0.5 * C_Y.col(nu).transpose() * dt_mus[mu.bfoff_in_shell].col(Y.bfoff_in_shell).segment(
1370  // obs->shell_to_function(obs->shell_on_center(Y.center, 0)),
1371  // obs->nbasis_on_center(Y.center)
1372  // );
1373  // }
1374  //}
1375  */
1376  //----------------------------------------//
1377  }
1378  if(ish != jsh){
1379  for(auto rho : function_range(jsh)) {
1380  // dt_rhos[rho.bfoff_in_shell] is (nbf x Ysh.nbf)
1381  // C_Y is (Y.{obs_}atom_nbf x nbf)
1382  // result should be (Y.{obs_}atom_nbf x 1)
1383  Kt_part.row(rho).segment(obs_atom_bfoff, obs_atom_nbf).transpose() -=
1384  0.5 * C_Y * dt_rhos[rho.bfoff_in_shell].col(Y.bfoff_in_block);
1385  // The sigma <-> nu term
1386  Kt_part.row(rho).transpose() -= 0.5
1387  * C_Y.transpose()
1388  * dt_rhos[rho.bfoff_in_shell].col(Y.bfoff_in_block).segment(obs_atom_bfoff, obs_atom_nbf);
1389  // Add back in the nu.center == Y.center part
1390  Kt_part.row(rho).segment(obs_atom_bfoff, obs_atom_nbf).transpose() += 0.5
1391  * C_Y.middleCols(obs_atom_bfoff, obs_atom_nbf).transpose()
1392  * dt_rhos[rho.bfoff_in_shell].col(Y.bfoff_in_block).segment(obs_atom_bfoff, obs_atom_nbf);
1393  //----------------------------------------//
1394  /* Failsafe version:
1395  //for(auto nu : function_range(obs)) {
1396  // for(auto sigma : function_range(obs)) {
1397  // if(nu.center == Y.center){
1398  // Kt_part(rho, nu) -= 0.5 * C_Y(nu.bfoff_in_atom, sigma) * dt_rhos[rho.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1399  // }
1400  // else if(sigma.center == Y.center){
1401  // Kt_part(rho, nu) -= 0.5 * C_Y(sigma.bfoff_in_atom, nu) * dt_rhos[rho.bfoff_in_shell](sigma, Y.bfoff_in_shell);
1402  // }
1403  // }
1404  //}
1405  */
1406  //----------------------------------------//
1407  }
1408  }
1409  }
1410  //----------------------------------------//
1411  /* Old version
1412  for(auto ksh : shell_range(obs, dfbs_)) {
1413  for(auto lsh : shell_range(obs, dfbs_, 0, ksh)) {
1414  //for(auto lsh : iter_significant_partners(ksh)) {
1415  if(ksh.center != Ysh.center and lsh.center != Ysh.center) continue;
1416  //----------------------------------------//
1417  //const double pf = ksh == lsh ? 1.0 : 2.0;
1418  //----------------------------------------//
1419  //for(int nu = ksh.bfoff; nu <= ksh.last_function; ++nu){ // function_range(ksh)) {
1420  // for(int sigma = lsh.bfoff; sigma <= lsh.last_function; ++sigma){ // function_range(ksh)) {
1421  for(auto nu : function_range(ksh)) {
1422  for(auto sigma : function_range(lsh)) {
1423  //----------------------------------------//
1424  // Get the coefficients
1425  IntPair nu_sigma(nu, sigma);
1426  auto cpair = coefs_[nu_sigma];
1427  VectorMap& C = (Ysh.center == ksh.center) ? *(cpair.first) : *(cpair.second);
1428  //----------------------------------------//
1429  //for(auto Y : function_range(Ysh)) {
1430  // auto& C_Y = coefs_transpose_[Y];
1431  // for(auto mu : function_range(ish)) {
1432  // Kt_part(mu, nu) -= 0.5
1433  // * dt_mus[mu.bfoff_in_shell].row(sigma)
1434  // * C_Y(nu.bfoff_in_atom, sigma);
1435  // }
1436  //}
1437  //----------------------------------------//
1438  for(auto mu : function_range(ish)) {
1439  Kt_part(mu, nu) -= 0.5
1440  * dt_mus[mu.bfoff_in_shell].row(sigma)
1441  * C.segment(Ysh.bfoff_in_atom, Ysh.nbf);
1442  if(ksh != lsh){
1443  Kt_part(mu, sigma) -= 0.5
1444  * dt_mus[mu.bfoff_in_shell].row(nu)
1445  * C.segment(Ysh.bfoff_in_atom, Ysh.nbf);
1446  }
1447  } // end loop over mu in ish
1448  //----------------------------------------//
1449  if(ish != jsh){
1450  for(auto rho : function_range(jsh)) {
1451  Kt_part(rho, nu) -= 0.5
1452  * dt_rhos[rho.bfoff_in_shell].row(sigma)
1453  * C.segment(Ysh.bfoff_in_atom, Ysh.nbf);
1454  if(ksh != lsh){
1455  Kt_part(rho, sigma) -= 0.5
1456  * dt_rhos[rho.bfoff_in_shell].row(nu)
1457  * C.segment(Ysh.bfoff_in_atom, Ysh.nbf);
1458  }
1459  } // end loop over rho in jsh
1460  } // end if ish != jsh
1461  //----------------------------------------//
1462  } // end loop over sigma in lsh
1463  } // end loop over nu in ksh
1464  //----------------------------------------//
1465  } // end loop over lsh
1466  } // end loop over ksh
1467  */
1468  //----------------------------------------//
1469  if(ithr == 0) timer.exit();
1470  /********************************************************/ #endif //2}}}
1471  /*-----------------------------------------------------*/
1472  } // end loop over shells Ysh
1473  } // end while get_shell_pair(ish, jsh)
1474  #endif
1475 
1477 // other stuff
1478 
1479  /*
1480  shell_iter_arbitrary_wrapper<std::vector<int>>
1481  iter_significant_partners(
1482  const ShellData& ish
1483  )
1484  {
1485  return shell_iter_arbitrary_wrapper<std::vector<int>>(
1486  shell_to_sig_shells_[ish.index],
1487  ish.basis,
1488  ish.dfbasis
1489  );
1490  }
1491  */
1492 
1493 
1494 // From J
1495  /*
1496  for(auto mu : function_range(ish)){
1497  for(auto nu : function_range(jsh)){
1498  //for(int jbf = 0, nu = jsh.bfoff; jbf < jsh.nbf; ++jbf, ++nu){
1499  //----------------------------------------//
1500  const int ijbf = mu.bfoff_in_shell*jsh.nbf + nu.bfoff_in_shell;
1501  //----------------------------------------//
1502  // compute dtilde contribution
1503  //dt.segment(Xblk.bfoff, Xblk.nbf) += perm_fact * D(mu, nu) * g_part->row(ijbf);
1504  //----------------------------------------//
1505  // add C_tilde * g_part to J
1506  if(nu <= mu){
1507  // only constructing the lower triangle of the J matrix
1508  jpart(mu, nu) += g_part->row(ijbf) * C_tilde.segment(Xblk.bfoff, Xblk.nbf);
1509  }
1510  //----------------------------------------//
1511  } // end loop over jbf
1512  } // end loop over ibf
1513  */
1514  //for(auto kshdf : shell_range(dfbs_)){
1515  // std::shared_ptr<Eigen::MatrixXd> g_part = ints_to_eigen(
1516  // ish, jsh, kshdf,
1517  // eris_3c_[ithr],
1518  // coulomb_oper_type_
1519  // );
1520  // for(int ibf = 0, mu = ish.bfoff; ibf < ish.nbf; ++ibf, ++mu){
1521  // //for(auto jbf : function_range(jsh)){
1522  // for(int jbf = 0, nu = jsh.bfoff; jbf < jsh.nbf; ++jbf, ++nu){
1523  // //----------------------------------------//
1524  // const int ijbf = ibf*jsh.nbf + jbf;
1525  // //----------------------------------------//
1526  // // compute dtilde contribution
1527  // dt.segment(kshdf.bfoff, kshdf.nbf) += perm_fact * D(mu, nu) * g_part->row(ijbf);
1528  // //----------------------------------------//
1529  // // add C_tilde * g_part to J
1530  // if(nu <= mu){
1531  // // only constructing the lower triangle of the J matrix
1532  // jpart(mu, nu) += g_part->row(ijbf) * C_tilde.segment(kshdf.bfoff, kshdf.nbf);
1533  // }
1534  // //----------------------------------------//
1535  // } // end loop over jbf
1536  // } // end loop over ibf
1537  //} // end loop over kshbf
1538 
1539 #endif /* CADF_ATTIC_H_ */
1540 
ShellBlockIterator< Iterable >
detail::basis_iterable
Definition: cadf_attic.h:222
mpqc::ci::sigma
void sigma(const CI< Type, Index > &ci, const mpqc::Vector &h, const Matrix &V, ci::Vector &C, ci::Vector &S)
Computes sigma 1,2,3 contributions.
Definition: sigma.hpp:30
shell_block_iterable
Definition: cadf_attic.h:406
ShellData
Definition: cadf_attic.h:100
ShellBlockIterator::Skeleton
Definition: cadf_attic.h:595
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
IterableBasisElementData
Definition: cadf_attic.h:445
function_iterable
Definition: cadf_attic.h:375
IterableBasisElementData< Iterator, ShellBlockIterator< Iterable > >
Definition: cadf_attic.h:630

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