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) \
43 inline const type BOOST_PP_CAT(get_, name)() const { assert(condition); return getter }\
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)\
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)\
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)
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) }
67 #define PROP_INIT2(r,data,impl) BOOST_PP_SEQ_ELEM(0, impl) = BOOST_PP_SEQ_ELEM(1, impl)
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)\
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) \
83 BOOST_PP_SEQ_FOR_EACH_R(1, DEF_PROPERTY, (clname)(type)(true), props) \
85 BOOST_PP_IF(BOOST_PP_CAT(clname, _USE_PROPERTIES), \
87 void init_function_name() { \
88 BOOST_PP_SEQ_FOR_EACH_I(PROPS_ASSIGN_impl, clname, props) \
93 void init_function_name() { \
94 BOOST_PP_SEQ_FOR_EACH_I( PROP_INIT, props, props) \
99 #define ShellData_USE_PROPERTIES 1
105 basis->shell_to_function(index);
108 basis->shell(index).nfunction();
111 basis->shell_to_center(index);
114 basis->shell_on_center(center, 0);
117 basis->shell_to_function(atom_shoff);
120 basis->nshell_on_center(center);
123 basis->nbasis_on_center(center);
131 ((atom_last_function)(
132 atom_bfoff + atom_nbf - 1;
135 atom_shoff + atom_nsh - 1;
140 (CONDITIONAL_PROPERTIES(dfbasis != 0,
142 dfbasis->shell_on_center(center, 0);
145 dfbasis->shell_to_function(atom_dfshoff);
148 dfbasis->nshell_on_center(center);
151 dfbasis->nbasis_on_center(center);
153 ((atom_df_last_function)(
154 atom_dfbfoff + atom_dfnbf - 1; ))
155 ((atom_df_last_shell)(
156 atom_dfshoff + atom_dfnsh - 1;
163 : BasisElementData(NotAssigned, 0, 0)
168 GaussianBasisSet* basis,
169 GaussianBasisSet* dfbasis = 0
170 ) : BasisElementData(idx, basis, dfbasis)
185 #if !ShellData_USE_PROPERTIES
191 void set_index(
int idx) {
193 #if !ShellData_USE_PROPERTIES
198 const ShellData& operator*()
const {
return *
this; }
205 dfbasis =
other.dfbasis;
212 operator int() { ASSERT_SHELL_BOUNDS;
return index; }
213 operator const int()
const { ASSERT_SHELL_BOUNDS;
return index; }
221 template <
typename DataContainer,
typename Iterable>
226 enum { NoLastIndex = -1, NoFirstIndex = -2 };
228 GaussianBasisSet* basis_;
229 GaussianBasisSet* dfbasis_;
230 const Iterable& indices_;
234 std::shared_ptr<Iterable> created_index_iterator_;
239 typedef decltype(Iterable().begin()) index_iterator_type;
245 GaussianBasisSet* basis,
246 GaussianBasisSet* dfbasis,
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
255 indices_(*created_index_iterator_)
259 const Ref<GaussianBasisSet>& basis,
260 const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
264 template<
typename iterable_like,
typename=
typename std::enable_if<
265 std::is_base_of<Iterable, typename std::decay<iterable_like>::type>::value
269 const Ref<GaussianBasisSet>& basis,
270 const OptionalRefParameter<GaussianBasisSet>& dfbasis,
271 const iterable_like& indices
272 ) : basis_(basis), dfbasis_(dfbasis), indices_(indices)
275 template<
typename iterable_like,
typename=
typename std::enable_if<
276 std::is_base_of<Iterable, typename std::decay<iterable_like>::type>::value
280 const Ref<GaussianBasisSet>& basis,
281 const iterable_like& indices,
282 const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
283 ) : basis_(basis), dfbasis_(dfbasis), indices_(indices)
287 const Ref<GaussianBasisSet>& basis,
289 const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
293 template<
typename int_like,
typename=
typename std::enable_if<
294 std::is_convertible<int_like, int>::value
298 const Ref<GaussianBasisSet>& basis,
299 const Ref<GaussianBasisSet>& dfbasis,
305 const Ref<GaussianBasisSet>& basis,
308 const OptionalRefParameter<GaussianBasisSet>& dfbasis = 0
313 const DataContainer& first,
314 const DataContainer& last
315 ) :
basis_iterable(first.basis, first.dfbasis, Iterable(first, last))
329 template <
typename Iterable>
336 typedef decltype(Iterable().begin()) index_iterator_type;
353 template<
typename ShellIncrementable>
355 make_shell_iterable(ShellIncrementable first, ShellIncrementable last) {
357 static_assert(not std::is_convertible<ShellIncrementable, BasisFunctionData>::value,
358 "Can't construct a shell_iterable using BasisFunctionData objects."
364 boost::iterator_range<boost::counting_iterator<ShellIncrementable>>(
374 template <
typename Iterable =
int_range>
377 static constexpr
int NotAssigned = -1;
379 int block_offset = NotAssigned;
385 typedef decltype(Iterable().begin()) index_iterator_type;
405 template<
typename Iterable>
408 int target_size = DEFAULT_TARGET_BLOCK_SIZE;
411 int reqs = SameCenter;
417 typedef decltype(Iterable().begin()) index_iterator_type;
432 target_size = new_target_size;
444 template<
typename Iterator,
typename ParentContainer>
452 typedef boost::forward_traversal_tag iterator_category;
453 typedef typename Iterator::difference_type difference_type;
455 template<
typename... Args>
458 GaussianBasisSet* basis,
459 GaussianBasisSet* dfbasis,
461 ) : spot(index_iter), ParentContainer(*index_iter, basis, dfbasis, std::forward<Args>(args)...)
467 ParentContainer::index = *spot;
468 ParentContainer::init();
478 template<
typename range_type>
483 static constexpr
int NoMaximumBlockSize = -20;
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()),
498 shell_iter(first_shell, last_shell)
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),
513 first_shell(possible_shell_iter.begin()),
514 last_shell(possible_shell_iter.end()),
515 shell_iter(first_shell, last_shell)
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
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),
534 first_shell(possible_shell_iter.begin()),
535 last_shell(possible_shell_iter.end()),
536 shell_iter(first_shell, last_shell)
551 return first_shell.index !=
other.first_shell.index or nshell !=
other.nshell;
556 GaussianBasisSet* basis;
557 GaussianBasisSet* dfbasis;
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;
579 int atom_df_last_function = -1;
580 int atom_df_last_shell = -1;
582 static int max_index(GaussianBasisSet* basis) {
return basis->nshell(); }
602 GaussianBasisSet* basis;
603 GaussianBasisSet* dfbasis;
607 ) : first_index(data.first_index),
610 dfbasis(data.dfbasis),
615 return first_index <
other.first_index
616 or (first_index ==
other.first_index and size <
other.size);
620 return first_index !=
other.first_index or size !=
other.size;
629 template<
typename Iterator,
typename Iterable>
654 template<
typename... Args>
656 const Iterable& index_iter,
657 GaussianBasisSet* basis,
658 GaussianBasisSet* dfbasis,
660 ) : spot(index_iter.begin()),
664 std::forward<Args>(args)...
672 const self_type& advance_to_last_block();
676 template <
typename DataContainer,
typename Iterable>
680 return value_type(indices_.begin(), basis_, dfbasis_);
683 template <
typename DataContainer,
typename Iterable>
687 auto last = indices_.end();
689 return value_type(last, basis_, dfbasis_);
692 template <
typename DataContainer,
typename Iterable>
696 return value_type(indices_.end(), basis_, dfbasis_);
699 template <
typename Iterable>
705 base_type::basis_, base_type::dfbasis_,
710 template <
typename Iterable>
715 Iterable(base_type::indices_.end(), base_type::indices_.end()),
716 base_type::basis_, base_type::dfbasis_,
721 template <
typename Iterable>
726 base_type::indices_.begin(),
727 base_type::indices_.end(),
728 base_type::basis_, base_type::dfbasis_,
731 rv.advance_to_last_block();
735 template<
typename Iterator,
typename Iterable>
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.",
748 auto iter_tmp = possible_shell_iter;
749 const auto& last_possible_shell = possible_shell_iter.last();
751 while(iter_tmp.begin() != iter_tmp.end()) {
752 first_shell = iter_tmp.begin();
756 const int first_center = first_shell.center;
757 auto& first_am = basis->shell(first_shell).am();
758 for(
auto ish : possible_shell_iter){
761 ((reqs & SameCenter) and ish.center != first_center)
764 ((reqs & SameAngularMomentum) and
765 basis->shell(ish).am() != first_am) or
767 (target_size != NoMaximumBlockSize and nbf >= target_size)
783 auto new_first_shell = last_shell;
785 iter_tmp = make_shell_iterable(new_first_shell, last_possible_shell);
788 possible_shell_iter = make_shell_iterable(first_shell, last_possible_shell);
792 template<
typename Iterable>
796 ) :
ShellBlockIterator<Iterable>(sk.first_index, sk.basis, sk.dfbasis, sk.reqs, sk.size)
800 template<
typename Iterator,
typename Iterable>
805 possible_shell_iter = make_shell_iter(spot, possible_shell_iter.last());
858 for(
auto lsh : iter_significant_partners(ksh)) {
859 if(ksh.center != Xblk.center && lsh.center != Xblk.center)
continue;
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;
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;
878 Kt_part(mu, nu) += C[X.bfoff_in_atom] * B_mus[mu.bfoff_in_shell](
sigma, X.bfoff_in_block);
890 while(get_shell_pair(ish, jsh, SignificantPairs)){
891 const double pf = 2.0;
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);
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);
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);
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();
933 if(ithr == 0) timer.enter(
"compute ints");
934 std::shared_ptr<Eigen::MatrixXd> g_part = ints_to_eigen(
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);
948 if(ithr == 0) timer.change(
"B_rhos");
951 for(
int ibf = 0, mu = ish.bfoff; ibf < ish.nbf; ++ibf, ++mu){
952 for(
int jbf = 0; jbf < jsh.nbf; ++jbf){
954 B_rhos[jbf].middleCols(kshdf_bfoff, kshdf_nbf) +=
955 pf * D.col(mu) * g_part->row(ibf * jsh.nbf + jbf);
961 if(ithr == 0) timer.exit();
966 write_as_xml(
"B_part", B_mus[mu.bfoff_in_shell], std::map<std::string, int>{ {
"mu", mu}, {
"jsh", jsh} } );
970 write_as_xml(
"B_part", B_rhos[rho.bfoff_in_shell], std::map<std::string, int>{ {
"mu", rho}, {
"jsh", jsh} } );
975 if(ithr == 0) timer.exit(
"build B");
1077 if(ithr == 0) timer.enter(
"K contributions");
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);
1087 Kt_part.row(mu).segment(obs_atom_bfoff, obs_atom_nbf).transpose() +=
1088 C_Y * B_mus[mu.bfoff_in_shell].col(Y);
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);
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);
1124 Kt_part.row(rho).segment(obs_atom_bfoff, obs_atom_nbf).transpose() +=
1125 C_Y * B_rhos[rho.bfoff_in_shell].col(Y);
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);
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);
1201 if(ithr == 0) timer.exit(
"K contributions");
1231 while(get_shell_pair(ish, jsh, SignificantPairs)){
1234 std::vector<Eigen::MatrixXd> Ct_mus(ish.nbf);
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);
1241 if(ithr == 0) timer.enter(
"form Ct_mu");
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(
1250 if(ithr == 0) timer.exit(
"compute ints");
1251 const auto& g2_part = *g2_part_ptr;
1256 IntPair mu_rho(mu, rho);
1257 auto cpair = coefs_[mu_rho];
1258 VectorMap& Ca = *(cpair.first);
1260 Ct_mus[mu.bfoff_in_shell].row(rho.bfoff_in_shell).transpose() +=
1261 g2_part * Ca.segment(Xblk.bfoff_in_atom, Xblk.nbf);
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(
1274 const auto& g2_part = *g2_part_ptr;
1279 IntPair mu_rho(mu, rho);
1280 auto cpair = coefs_[mu_rho];
1281 VectorMap& Cb = *(cpair.second);
1283 Ct_mus[mu.bfoff_in_shell].row(rho.bfoff_in_shell).transpose() +=
1284 g2_part * Cb.segment(Xblk.bfoff_in_atom, Xblk.nbf);
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);
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];
1314 dt_rhos[rho.bfoff_in_shell].resize(nbf, Yblk.nbf);
1315 dt_rhos[rho.bfoff_in_shell] = Eigen::MatrixXd::Zero(nbf, Yblk.nbf);
1317 dt_rhos[rho.bfoff_in_shell] += 2.0 * D.col(mu) * Ct_mus[mu.bfoff_in_shell].row(rho.bfoff_in_shell);
1335 if(ithr == 0) timer.change(
"K contributions");
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);
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);
1348 Kt_part.row(mu).transpose() -= 0.5
1350 * dt_mus[mu.bfoff_in_shell].col(Y.bfoff_in_block).segment(obs_atom_bfoff, obs_atom_nbf);
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);
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);
1386 Kt_part.row(rho).transpose() -= 0.5
1388 * dt_rhos[rho.bfoff_in_shell].col(Y.bfoff_in_block).segment(obs_atom_bfoff, obs_atom_nbf);
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);
1469 if(ithr == 0) timer.exit();