MPQC  3.0.0-alpha
fockbuilder.h
1 //
2 // fockbuilder.h
3 //
4 // Copyright (C) 2009 Edward Valeev
5 //
6 // Author: Edward Valeev <evaleev@vt.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifndef _mpqc_src_lib_chemistry_qc_lcao_fockbuilder_h
29 #define _mpqc_src_lib_chemistry_qc_lcao_fockbuilder_h
30 
31 #include <cassert>
32 #include <util/ref/ref.h>
33 #include <util/group/thread.h>
34 #include <util/misc/scexception.h>
35 #include <chemistry/qc/basis/gpetite.h>
36 #include <chemistry/qc/basis/symmint.h>
37 #include <chemistry/qc/lcao/fockbuild.h>
38 #include <chemistry/qc/lcao/clhfcontrib.h>
39 #include <chemistry/qc/lcao/hsoshfcontrib.h>
40 #include <chemistry/qc/lcao/moints_runtime.h>
41 #include <chemistry/qc/lcao/fockbuild_runtime.h>
42 
43 namespace sc {
44 
45  namespace detail {
46  template<bool bra_eq_ket> struct FockMatrixType;
47  template<> struct FockMatrixType<true> {
48  typedef RefSymmSCMatrix value;
49  struct Factory {
50  static value create(const Ref<SCMatrixKit>& kit, const RefSCDimension& bradim, const RefSCDimension& ketdim) {
51  return kit->symmmatrix(bradim);
52  }
53  static void transform(const value& result, const value& original,
54  const RefSCMatrix& Ubra, const RefSCMatrix& Uket,
55  SCMatrix::Transform transform_type = SCMatrix::NormalTransform) {
56  result.assign(0.0);
57  result.accumulate_transform(Ubra,original,transform_type);
58  }
59  static void convert(const value& result, const value& original) {
60  return result->convert(original);
61  }
62  template <typename T> static void convert(const value& result, const T& original) {
63  abort();
64  }
65  };
66  };
67  template<> struct FockMatrixType<false> {
68  typedef RefSCMatrix value;
69  struct Factory {
70  static value create(const Ref<SCMatrixKit>& kit, const RefSCDimension& bradim, const RefSCDimension& ketdim) {
71  return kit->matrix(bradim,ketdim);
72  }
73  static void transform(const value& result, const value& original,
74  const RefSCMatrix& Ubra, const RefSCMatrix& Uket,
75  SCMatrix::Transform transform_type = SCMatrix::NormalTransform) {
76  if (transform_type == SCMatrix::NormalTransform)
77  result.accumulate_product(Ubra, original * Uket.t());
78  else
79  result.accumulate_product(Ubra.t(), original * Uket);
80  }
81  static void convert(const value& result, const value& original) {
82  return result->convert(original);
83  }
84  template <typename T> static void convert(const value& result, const T& original) {
85  abort();
86  }
87  };
88  };
89 
90  typedef Ref<OneBodyInt> (Integral::*OneBodyIntCreator)();
91  typedef Ref<OneBodyInt> (Integral::*MultipoleOneBodyIntCreator)(const Ref<IntParamsOrigin>&);
92 
93  template <OneBodyIntCreator OpEval>
94  RefSymmSCMatrix onebodyint(const Ref<GaussianBasisSet>& bas,
95  const Ref<Integral>& integral) {
96  Ref<Integral> localints = integral->clone();
97  localints->set_basis(bas);
98  Ref<PetiteList> pl = localints->petite_list();
99 
100  // form skeleton operator matrix in AO basis
101  RefSymmSCMatrix opmat_ao(bas->basisdim(), bas->matrixkit());
102  opmat_ao.assign(0.0);
103  Ref<SCElementOp> elemop =
104  new OneBodyIntOp(new SymmOneBodyIntIter( (localints->*OpEval)(), pl));
105  opmat_ao.element_op(elemop);
106  elemop= 0;
107 
108  // now symmetrize
109  RefSymmSCMatrix opmat(pl->SO_basisdim(), bas->so_matrixkit());
110  pl->symmetrize(opmat_ao, opmat);
111 
112  return opmat;
113  }
114 
115  template <OneBodyIntCreator OpEval>
116  RefSymmSCMatrix onebodyint_ao(const Ref<GaussianBasisSet>& bas,
117  const Ref<Integral>& integral) {
118  RefSymmSCMatrix opmat_so = onebodyint<OpEval>(bas, integral);
119  Ref<Integral> localints = integral->clone();
120  localints->set_basis(bas);
121  Ref<PetiteList> pl = localints->petite_list();
122 
123  RefSymmSCMatrix opmat_ao = pl->to_AO_basis(opmat_so);
124  return opmat_ao;
125  }
126 
127  template <OneBodyIntCreator OpEval>
128  RefSCMatrix onebodyint_ao(const Ref<GaussianBasisSet>& brabas,
129  const Ref<GaussianBasisSet>& ketbas,
130  const Ref<Integral>& integral) {
131 
132  Ref<Integral> localints = integral->clone();
133  Ref<GPetiteList2> pl12 = GPetiteListFactory::plist2(brabas,ketbas);
134  localints->set_basis(brabas,ketbas);
135  Ref<OneBodyInt> ov_ints = (localints->*OpEval)();
136 
137  // form overlap in AO basis
138  RefSCMatrix sao(brabas->basisdim(), ketbas->basisdim(), brabas->matrixkit());
139  sao.assign(0.0);
140 
141  // TODO make more efficient and parallelizable by switching to operators
142  const Ref<GaussianBasisSet>& bs1 = brabas;
143  const Ref<GaussianBasisSet>& bs2 = ketbas;
144  const int nshell1 = bs1->nshell();
145  const int nshell2 = bs2->nshell();
146  for(int sh1=0; sh1<nshell1; sh1++) {
147  const int bf1_offset = bs1->shell_to_function(sh1);
148  const int nbf1 = bs1->shell(sh1).nfunction();
149 
150  for(int sh2=0; sh2<nshell2; sh2++) {
151  const int bf2_offset = bs2->shell_to_function(sh2);
152  const int nbf2 = bs2->shell(sh2).nfunction();
153 
154  ov_ints->compute_shell(sh1,sh2);
155  const double *ovintsptr = ov_ints->buffer();
156 
157  int bf1_index = bf1_offset;
158  for(int bf1=0; bf1<nbf1; bf1++, bf1_index++, ovintsptr+=nbf2) {
159  int bf2_index = bf2_offset;
160  const double *ptr = ovintsptr;
161  for(int bf2=0; bf2<nbf2; bf2++, bf2_index++) {
162 
163  sao.set_element(bf1_index, bf2_index, *(ptr++));
164 
165  }
166  }
167  }
168  }
169 
170  return sao;
171 
172  }
173 
174  template <MultipoleOneBodyIntCreator OpEval>
175  void onebodyint_ao(const Ref<GaussianBasisSet>& brabas,
176  const Ref<GaussianBasisSet>& ketbas,
177  const Ref<Integral>& integral, const Ref<IntParamsOrigin>& multipole_origin,
178  std::vector<RefSCMatrix>& result) {
179 
180  Ref<Integral> localints = integral->clone();
181  Ref<GPetiteList2> pl12 = GPetiteListFactory::plist2(brabas,ketbas);
182  localints->set_basis(brabas,ketbas);
183  Ref<OneBodyInt> ob_ints = (localints->*OpEval)(multipole_origin);
184 
185  // form obints in AO basis
186  const size_t nopers = result.size();
187  for(int i=0; i<nopers; ++i) {
188  result[i] = brabas->matrixkit()->matrix(brabas->basisdim(), ketbas->basisdim());
189  result[i].assign(0.0);
190  }
191 
192  const Ref<GaussianBasisSet>& bs1 = brabas;
193  const Ref<GaussianBasisSet>& bs2 = ketbas;
194  const int nshell1 = bs1->nshell();
195  const int nshell2 = bs2->nshell();
196  for(int sh1=0; sh1<nshell1; sh1++) {
197  const int bf1_offset = bs1->shell_to_function(sh1);
198  const int nbf1 = bs1->shell(sh1).nfunction();
199 
200  for(int sh2=0; sh2<nshell2; sh2++) {
201  const int bf2_offset = bs2->shell_to_function(sh2);
202  const int nbf2 = bs2->shell(sh2).nfunction();
203 
204  ob_ints->compute_shell(sh1,sh2);
205  const double *ovintsptr = ob_ints->buffer();
206 
207  int bf1_index = bf1_offset;
208  for(int bf1=0; bf1<nbf1; bf1++, bf1_index++, ovintsptr+=nopers*nbf2) {
209  int bf2_index = bf2_offset;
210  const double *ptr = ovintsptr;
211  for(int bf2=0; bf2<nbf2; bf2++, bf2_index++) {
212 
213  for(size_t oper=0; oper<nopers; ++oper)
214  result[oper].set_element(bf1_index, bf2_index, *(ptr++));
215 
216  }
217  }
218  }
219  }
220 
221  }
222 
223  template <OneBodyIntCreator OpEval>
224  RefSCMatrix onebodyint(const Ref<GaussianBasisSet>& brabas,
225  const Ref<GaussianBasisSet>& ketbas,
226  const Ref<Integral>& integral) {
227 
228  RefSCMatrix sao = onebodyint_ao<OpEval>(brabas, ketbas, integral);
229  Ref<Integral> braints = integral->clone(); braints->set_basis(brabas);
230  Ref<PetiteList> brapl = braints->petite_list();
231  Ref<Integral> ketints = integral->clone(); ketints->set_basis(ketbas);
232  Ref<PetiteList> ketpl = ketints->petite_list();
233 
234  // SO basis is always blocked, so first make sure a is blocked
235  RefSCMatrix sao_blk = dynamic_cast<BlockedSCMatrix*>(sao.pointer());
236  if (sao_blk.null()) {
237  sao_blk = brabas->so_matrixkit()->matrix(brapl->AO_basisdim(),ketpl->AO_basisdim());
238  sao_blk->convert(sao);
239  }
240  // if C1, then do nothing
241  if (brabas->molecule()->point_group()->order() == 1)
242  return sao_blk;
243 
244  // convert to SO basis
245  RefSCMatrix ssoao = brapl->aotoso().t() * sao_blk;
246  RefSCMatrix result = ssoao * ketpl->aotoso();
247 
248  return result;
249  }
250 
252  RefSymmSCMatrix overlap(const Ref<GaussianBasisSet>& bas,
253  const Ref<Integral>& integral);
255  RefSCMatrix overlap(const Ref<GaussianBasisSet>& brabs,
256  const Ref<GaussianBasisSet>& ketbs,
257  const Ref<Integral>& integral);
259  void edipole(const Ref<GaussianBasisSet>& basis,
260  const Ref<Integral>& integral,
261  RefSymmSCMatrix& mu_x,
262  RefSymmSCMatrix& mu_y,
263  RefSymmSCMatrix& mu_z
264  );
266  RefSymmSCMatrix twobody_twocenter_coulomb(const Ref<GaussianBasisSet>& bas,
267  const Ref<Integral>& integral);
269  RefSymmSCMatrix nonrelativistic(const Ref<GaussianBasisSet>& bas,
270  const Ref<Integral>& integral);
272  RefSymmSCMatrix pauli(const Ref<GaussianBasisSet>& bas,
273  const Ref<Integral>& integral);
275  RefSymmSCMatrix dk(int dklev, const Ref<GaussianBasisSet>& bs,
276  const Ref<GaussianBasisSet>& pbs,
277  const Ref<Integral>& integral);
278 
279  RefSCMatrix coulomb(const Ref<TwoBodyFourCenterMOIntsRuntime>& ints_rtime,
280  const Ref<OrbitalSpace>& occspace,
281  const Ref<OrbitalSpace>& braspace,
282  const Ref<OrbitalSpace>& ketspace);
283  RefSCMatrix exchange(const Ref<TwoBodyFourCenterMOIntsRuntime>& ints_rtime,
284  const Ref<OrbitalSpace>& occspace,
285  const Ref<OrbitalSpace>& braspace,
286  const Ref<OrbitalSpace>& ketspace);
287 
288  RefSCMatrix coulomb_df(const Ref<DensityFittingInfo>& df_info,
289  const RefSymmSCMatrix& P,
290  const Ref<GaussianBasisSet>& brabs,
291  const Ref<GaussianBasisSet>& ketbs,
292  const Ref<GaussianBasisSet>& obs);
293 #ifdef MPQC_NEW_FEATURES
294  RefSCMatrix coulomb_df_local(const Ref<DensityFittingInfo>& df_info,
295  const RefSymmSCMatrix& P,
296  const Ref<GaussianBasisSet>& brabs,
297  const Ref<GaussianBasisSet>& ketbs,
298  const Ref<GaussianBasisSet>& obs);
299 #endif // MPQC_NEW_FEATURES
300  RefSCMatrix exchange_df(const Ref<DensityFittingInfo>& df_info,
301  const RefSymmSCMatrix& P,
302  SpinCase1 spin,
303  const Ref<GaussianBasisSet>& brabs,
304  const Ref<GaussianBasisSet>& ketbs,
305  const Ref<GaussianBasisSet>& obs,
306  const Ref<FockBuildRuntime::PSqrtRegistry>& psqrtregistry);
307 #ifdef MPQC_NEW_FEATURES
308  RefSCMatrix exchange_df_local(const Ref<DensityFittingInfo>& df_info,
309  const RefSymmSCMatrix& P,
310  SpinCase1 spin,
311  const Ref<GaussianBasisSet>& brabs,
312  const Ref<GaussianBasisSet>& ketbs,
313  const Ref<GaussianBasisSet>& obs,
314  const Ref<FockBuildRuntime::PSqrtRegistry>& psqrtregistry);
315 #endif // MPQC_NEW_FEATURES
316  } // end of namespace detail
317 
318 
320  template<bool bra_eq_ket> class TwoBodyFockMatrixBuilder: public RefCount {
321 
322  public:
323 
324  typedef typename detail::FockMatrixType<bra_eq_ket>::value ResultType;
325  typedef typename detail::FockMatrixType<bra_eq_ket>::Factory ResultFactory;
326 
327  TwoBodyFockMatrixBuilder(bool compute_F,
328  bool compute_J,
329  bool compute_K,
330  const Ref<GaussianBasisSet>& brabasis,
331  const Ref<GaussianBasisSet>& ketbasis,
332  const Ref<GaussianBasisSet>& densitybasis,
333  const RefSymmSCMatrix& density,
334  const RefSymmSCMatrix& openshelldensity,
335  const Ref<Integral>& integral,
336  const Ref<MessageGrp>& msg,
337  const Ref<ThreadGrp>& thr,
338  double accuracy = 1e-20) :
339  compute_F_(compute_F),
340  compute_J_(compute_J),
341  compute_K_(compute_K)
342  {
343 
344  if (brabasis->equiv(ketbasis) != bra_eq_ket)
345  throw ProgrammingError("FockMatrixBuilder::FockMatrixBuilder -- inconsistent constructor and template arguments",
346  __FILE__, __LINE__);
347 
348  Ref<Integral> localints = integral->clone();
349 
350  Ref<FockContribution> fc;
351  const bool openshell = openshelldensity;
352  if (openshell) {
353  fc = new HSOSHFContribution(brabasis, ketbasis, densitybasis, std::string("replicated"));
354  ntypes_ = 2;
355  fc->set_pmat(0, density);
356  fc->set_pmat(1, openshelldensity);
357  } else {
358  fc = new CLHFContribution(brabasis, ketbasis, densitybasis, std::string("replicated"));
359  ntypes_ = 1;
360  fc->set_pmat(0, density);
361  }
362 
363  // FockBuild can compute either J and K separately, or the total F.
364  // Determine whether we really need to compute J and K, or F
365  const bool really_compute_F = compute_F_ && !compute_J_ && !compute_K_;
366  const bool really_compute_J = !really_compute_F && compute_J_;
367  const bool really_compute_K = !really_compute_F && compute_K_;
368 
369  // FockBuild only accepts matrices created with appropriate basisdim(), which are not blocked, hence must
370  // use non-blocked kit
371  ResultType G_ao_skel[2][3];
372  for(int t=0; t<ntypes_; ++t) {
373  if (really_compute_J) {
374  G_ao_skel[t][0] = ResultFactory::create(brabasis->matrixkit(),brabasis->basisdim(),ketbasis->basisdim());
375  G_ao_skel[t][0].assign(0.0);
376  fc->set_jmat(t, G_ao_skel[t][0]);
377  }
378  if (really_compute_K) {
379  G_ao_skel[t][1] = ResultFactory::create(brabasis->matrixkit(),brabasis->basisdim(),ketbasis->basisdim());
380  G_ao_skel[t][1].assign(0.0);
381  fc->set_kmat(t, G_ao_skel[t][1]);
382  }
383  if (really_compute_F) {
384  G_ao_skel[t][2] = ResultFactory::create(brabasis->matrixkit(),brabasis->basisdim(),ketbasis->basisdim());
385  G_ao_skel[t][2].assign(0.0);
386  fc->set_fmat(t, G_ao_skel[t][2]);
387  }
388  }
389 
390  const bool prefetch_blocks = false;
391  Ref<FockDistribution> fd = new FockDistribution;
392  fb_ = new FockBuild(fd, fc, prefetch_blocks, brabasis, ketbasis, densitybasis, msg, thr, localints);
393  fb_->set_accuracy(accuracy);
394  fb_->set_compute_J(really_compute_J);
395  fb_->set_compute_K(really_compute_K);
396  fb_->build();
397 
398  Ref<GPetiteList2> pl = GPetiteListFactory::plist2(brabasis,ketbasis);
399  localints->set_basis(brabasis);
400  Ref<PetiteList> brapl = localints->petite_list();
401  localints->set_basis(ketbasis);
402  Ref<PetiteList> ketpl = localints->petite_list();
403  const int ng = pl->point_group()->char_table().order();
404 
406  for(int t=0; t<ntypes_; ++t) {
407  for(int c=0; c<3; ++c) {
408 
409  if (really_compute_J == false && c == 0) continue;
410  if (really_compute_K == false && c == 1) continue;
411  if (really_compute_F == false && c == 2) continue;
412  // J matrices are spin-independent
413  if (c == 0 && t == 1) continue;
414 
415  // if C1 -- nothing else needs to be done, return the result
416  // same holds for brabasis != ketbasis since FockBuild does not use symmetry in that case
417  if (ng == 1 || !bra_eq_ket) {
418  RefSCDimension braaodim = brapl->AO_basisdim();
419  RefSCDimension ketaodim = ketpl->AO_basisdim();
420  result_[t][c] = ResultFactory::create(brabasis->so_matrixkit(),brapl->AO_basisdim(),ketpl->AO_basisdim());
421  result_[t][c]->convert(G_ao_skel[t][c]);
422  }
423  else { // if not C1 and , symmetrize the skeleton G matrix to produce the SO basis G matrix
424  G_ao_skel[t][c].scale(1.0/(double)ng);
425  ResultType G_so = ResultFactory::create(brabasis->so_matrixkit(),brapl->SO_basisdim(),ketpl->SO_basisdim());
426  symmetrize(pl,localints,G_ao_skel[t][c],G_so);
427  G_ao_skel[t][c] = 0;
428 
429  // and convert back to AO basis, but this time make a blocked matrix
430  RefSCDimension braaodim = brapl->AO_basisdim();
431  RefSCDimension ketaodim = ketpl->AO_basisdim();
432  result_[t][c] = ResultFactory::create(brabasis->so_matrixkit(),brapl->AO_basisdim(),ketpl->AO_basisdim());
433  ResultFactory::transform(result_[t][c], G_so, brapl->sotoao(), ketpl->sotoao(), SCMatrix::TransposeTransform);
434 
435  }
436 
437  // FockBuild computes -K, hence change the sign
438  if (c == 1) result_[t][1].scale(-1.0);
439 
440  }
441  }
442 
443  if (compute_F_ && !really_compute_F) {
444  // F = J - K = J - 1/2 (Ka + Kb)
445  result_[0][2] = result_[0][0] - result_[0][1];
446  if (ntypes_ == 2) { // Fo = -Ko = -1/2 (Ka - Kb)
447  result_[1][2] = result_[1][1].copy();
448  result_[1][2].scale(-1.0);
449  }
450  }
451 
452  }
453 
454  const Ref<FockBuild>& builder() const { return fb_; }
455  double nints() const { return builder()->contrib()->nint(); }
456  const ResultType& F(unsigned int t = 0) const {
457  MPQC_ASSERT(compute_F_ && t < ntypes_);
458  return result_[t][2];
459  }
460  const ResultType& J(unsigned int t = 0) const {
461  MPQC_ASSERT(compute_J_ && t < ntypes_ && t == 0);
462  return result_[t][0];
463  }
464  const ResultType& K(unsigned int t = 0) const {
465  MPQC_ASSERT(compute_K_ && t < ntypes_);
466  return result_[t][1];
467  }
468  ResultType F(SpinCase1 spin) const {
469  if (ntypes_ == 1)
470  return F(0);
471  else {
472  return (spin == Alpha) ? F(0) + F(1) : F(0) - F(1);
473  }
474  }
475  ResultType K(SpinCase1 spin) const {
476  if (ntypes_ == 1)
477  return K(0);
478  else {
479  return (spin == Alpha) ? K(0) + K(1) : K(0) - K(1);
480  }
481  }
482 
483  private:
484 
485  Ref<FockBuild> fb_;
486  bool compute_J_;
487  bool compute_K_;
488  bool compute_F_;
489  int ntypes_;
490  ResultType result_[2][3];
491 
492  }; // class TwoBodyFockMatrixBuilder
493 
495  class TwoBodyFockMatrixTransformBuilder: public RefCount {
496 
497  public:
498 
499  TwoBodyFockMatrixTransformBuilder(bool compute_J,
500  bool compute_K,
501  SpinCase1 spin,
502  const Ref<OrbitalSpace>& braspace,
503  const Ref<OrbitalSpace>& ketspace,
504  const Ref<OrbitalSpace>& occspace_A,
505  const Ref<OrbitalSpace>& occspace_B,
506  const Ref<TwoBodyFourCenterMOIntsRuntime>& ints_rtime) :
507  compute_J_(compute_J),
508  compute_K_(compute_K),
509  compute_F_(compute_J && compute_K)
510  {
511 
512  for(int t=0; t<3; ++t) {
513 
514  if (compute_J == false && t == 0) continue;
515  if (compute_K == false && t == 1) continue;
516 
517  if (t == 0) { // coulomb
518  result_[t] = detail::coulomb(ints_rtime,occspace_A,braspace,ketspace);
519  if (*occspace_A == *occspace_B) { // alpha and beta spins equivalent? Scale by 2, else compute
520  result_[t].scale(2.0);
521  }
522  else {
523  result_[t].accumulate( detail::coulomb(ints_rtime,occspace_B,braspace,ketspace) );
524  }
525  }
526 
527  if (t == 1) { // exchange
528  result_[t] = detail::exchange(ints_rtime,(spin == Alpha ? occspace_A : occspace_B),braspace,ketspace);
529  }
530 
531  }
532 
533  if (compute_F_)
534  result_[2] = result_[0] - result_[1];
535 
536  }
537 
538  const RefSCMatrix& F() const {
539  MPQC_ASSERT(compute_F_);
540  return result_[2];
541  }
542  const RefSCMatrix& J() const {
543  MPQC_ASSERT(compute_J_);
544  return result_[0];
545  }
546  const RefSCMatrix& K() const {
547  MPQC_ASSERT(compute_K_);
548  return result_[1];
549  }
550 
551  private:
552 
553  bool compute_J_;
554  bool compute_K_;
555  bool compute_F_;
556  RefSCMatrix result_[3];
557 
558  }; // class TwoBodyFockMatrixTransformBuilder
559 
561  class TwoBodyFockMatrixDFBuilder: public RefCount {
562 
563  typedef RefSCMatrix ResultType;
564 
565  public:
566 
567  TwoBodyFockMatrixDFBuilder(bool compute_F,
568  bool compute_J,
569  bool compute_K,
570  const Ref<GaussianBasisSet>& brabasis,
571  const Ref<GaussianBasisSet>& ketbasis,
572  const Ref<GaussianBasisSet>& densitybasis,
573  const RefSymmSCMatrix& density,
574  const RefSymmSCMatrix& openshelldensity,
575  const Ref<DensityFittingInfo>& df_info,
576  const Ref<FockBuildRuntime::PSqrtRegistry>& psqrtregistry) :
577  compute_J_(compute_J),
578  compute_K_(compute_K),
579  compute_F_(compute_F),
580  ntypes_(openshelldensity ? 2 : 1)
581  {
582 
583  // DF-based builds are separate for J and K separately
584  // Determine whether we really need to compute J and K, or F
585  const bool really_compute_F = compute_F_;
586  const bool really_compute_J = compute_J_ || compute_F_;
587  const bool really_compute_K = compute_K_ || compute_F_;
588 
589  for(int t=0; t<ntypes_; ++t) {
590  const SpinCase1 spin = static_cast<SpinCase1>(t);
591  for(int c=0; c<3; ++c) {
592 
593  if (really_compute_J == false && c == 0) continue;
594  if (really_compute_K == false && c == 1) continue;
595  if (really_compute_F == false && c == 2) continue;
596  // J matrices are spin-independent
597  if (c == 0 && t == 1) continue;
598 
599  if (c == 0) { // coulomb
600  if(df_info->params()->local_coulomb()) {
601 #ifdef MPQC_NEW_FEATURES
602  result_[0][c] = detail::coulomb_df_local(df_info, density, brabasis, ketbasis, densitybasis);
603 #else // MPQC_NEW_FEATURES
604  throw sc::FeatureNotImplemented("Local DF Coulomb build requires MPQC3 runtime: compile in boost, tiledarray, etc.",
605  __FILE__,__LINE__);
606 #endif // MPQC_NEW_FEATURES
607  }
608  else
609  result_[0][c] = detail::coulomb_df(df_info, density, brabasis, ketbasis, densitybasis);
610  }
611 
612  if (c == 1) { // exchange
613  RefSymmSCMatrix Pspin;
614  SpinCase1 spincase;
615  if (openshelldensity) {
616  Pspin = (spin == Alpha) ? density + openshelldensity : density - openshelldensity;
617  spincase = spin;
618  }
619  else {
620  Pspin = density.copy();
621  spincase = AnySpinCase1;
622  }
623  Pspin.scale(0.5);
624  if(df_info->params()->local_exchange()) {
625 #ifdef MPQC_NEW_FEATURES
626  result_[t][c] = detail::exchange_df_local(df_info, Pspin, spincase, brabasis, ketbasis, densitybasis,
627  psqrtregistry);
628 #else // MPQC_NEW_FEATURES
629  throw sc::FeatureNotImplemented("Local DF exchange build requires MPQC3 runtime: compile in boost, tiledarray, etc.",
630  __FILE__,__LINE__);
631 #endif // MPQC_NEW_FEATURES
632  }
633  else
634  result_[t][c] = detail::exchange_df(df_info, Pspin, spincase, brabasis, ketbasis, densitybasis,
635  psqrtregistry);
636  }
637 
638  if (c == 2) { // fock
639  result_[t][2] = result_[0][0] - result_[t][1];
640  }
641  }
642  }
643 
644  }
645 
646  const ResultType& F(unsigned int t = 0) const {
647  MPQC_ASSERT(compute_F_ && t < ntypes_);
648  return result_[t][2];
649  }
650  const ResultType& J(unsigned int t = 0) const {
651  MPQC_ASSERT(compute_J_ && t < ntypes_);
652  return result_[t][0];
653  }
654  const ResultType& K(unsigned int t = 0) const {
655  MPQC_ASSERT(compute_K_ && t < ntypes_);
656  return result_[t][1];
657  }
658  ResultType F(SpinCase1 spin) const {
659  if (ntypes_ == 1)
660  return F(0);
661  else {
662  return (spin == Alpha) ? F(0) : F(1);
663  }
664  }
665  ResultType K(SpinCase1 spin) const {
666  if (ntypes_ == 1)
667  return K(0);
668  else {
669  // DF-based computed alpha and beta exchange matrices
670  return (spin == Alpha) ? K(0) : K(1);
671  }
672  }
673 
674  private:
675 
676  bool compute_J_;
677  bool compute_K_;
678  bool compute_F_;
679  int ntypes_;
680  ResultType result_[2][3];
681 
682  }; // class TwoBodyFockMatrixDFBuilder
683 
685  template<bool bra_eq_ket> class OneBodyFockMatrixBuilder: public RefCount {
686 
687  public:
688 
689  typedef enum { NonRelativistic, Pauli, DK1, DK2 } OneBodyHamiltonianType;
690 
691  typedef typename detail::FockMatrixType<bra_eq_ket>::value ResultType;
692  typedef typename detail::FockMatrixType<bra_eq_ket>::Factory ResultFactory;
693 
694  OneBodyFockMatrixBuilder(OneBodyHamiltonianType type,
695  const Ref<GaussianBasisSet>& brabasis,
696  const Ref<GaussianBasisSet>& ketbasis,
697  const Ref<GaussianBasisSet>& pbasis,
698  const Ref<Integral>& integral,
699  double accuracy = 1e-12)
700  {
701  if (brabasis->equiv(ketbasis) != bra_eq_ket)
702  throw ProgrammingError("OneBodyFockMatrixBuilder::OneBodyFockMatrixBuilder -- inconsistent constructor and template arguments",
703  __FILE__, __LINE__);
704 
705  const Ref<GaussianBasisSet>& bs1 = brabasis;
706  const Ref<GaussianBasisSet>& bs2 = ketbasis;
707  const bool bs1_eq_bs2 = bra_eq_ket;
708 
709  Ref<GaussianBasisSet> hcore_basis = (bs1_eq_bs2) ? bs1 : bs1+bs2;
710 
711  Ref<GaussianBasisSet> p_basis = pbasis;
712  if (type == DK1 || type == DK2) {
713  // momentum basis in DKH calculations must span both bra and ket basis sets.
714  // easiest to achieve if both are included in p_basis
715  p_basis = p_basis + hcore_basis;
716  }
717 
718  RefSymmSCMatrix hsymm;
719  switch (type) {
720  case NonRelativistic: hsymm = detail::nonrelativistic(hcore_basis,integral); break;
721  case Pauli: hsymm = detail::pauli(hcore_basis,integral); break;
722 
723  int dklev;
724  case DK1: dklev = 1;
725  case DK2: dklev = 2;
726  hsymm = detail::dk(dklev, hcore_basis, p_basis, integral);
727  break;
728 
729  default:
730  throw ProgrammingError("Unrecognized Hamiltonian type",__FILE__,__LINE__);
731  }
732 
733  // convert hsymm to the AO basis
734  Ref<Integral> localints = integral->clone();
735  localints->set_basis(hcore_basis,hcore_basis);
736  Ref<PetiteList> hcore_pl = localints->petite_list();
737  RefSymmSCMatrix hsymm_ao = hcore_pl->to_AO_basis(hsymm);
738  hsymm = 0;
739 
740  localints->set_basis(brabasis);
741  Ref<PetiteList> brapl = localints->petite_list();
742  localints->set_basis(ketbasis);
743  Ref<PetiteList> ketpl = localints->petite_list();
744  RefSCDimension braaodim = brapl->AO_basisdim();
745  RefSCDimension ketaodim = ketpl->AO_basisdim();
746  if (bs1_eq_bs2) {
747  result_ = ResultFactory::create(brabasis->so_matrixkit(),braaodim,ketaodim);
748  ResultFactory::convert(result_,hsymm_ao);
749  }
750  else {
751  RefSCMatrix result_rect = brabasis->so_matrixkit()->matrix(braaodim,ketaodim);
752  RefSCMatrix hrect_ao = brabasis->so_matrixkit()->matrix(hsymm_ao.dim(),hsymm_ao.dim());
753  hrect_ao.assign(0.0);
754  hrect_ao.accumulate(hsymm_ao);
755 
756  // extract the bs1 by bs2 block:
757  // loop over all bs1 fblocks
758  // loop over all bs2 fblocks
759  // copy block to h
760  // end loop
761  // end loop
762  GaussianBasisSetMap bs1_to_hbs(bs1, hcore_basis);
763  GaussianBasisSetMap bs2_to_hbs(bs2, hcore_basis);
764  const int nfblock1 = bs1_to_hbs.nfblock();
765  const int nfblock2 = bs2_to_hbs.nfblock();
766  for (int rb = 0; rb < nfblock1; ++rb) {
767  const int rf_start1 = bs1_to_hbs.fblock_to_function(rb);
768  const int rf_end1 = rf_start1 + bs1_to_hbs.fblock_size(rb) - 1;
769 
770  const int rf_start12 = bs1_to_hbs.map_function(rf_start1);
771  const int rf_end12 = bs1_to_hbs.map_function(rf_end1);
772 
773  for (int cb = 0; cb < nfblock2; ++cb) {
774  const int cf_start2 = bs2_to_hbs.fblock_to_function(cb);
775  const int cf_end2 = cf_start2 + bs2_to_hbs.fblock_size(cb) - 1;
776 
777  const int cf_start12 =
778  bs2_to_hbs.map_function(cf_start2);
779  const int cf_end12 = bs2_to_hbs.map_function(cf_end2);
780 
781  // assign row-/col-subblock to h
782  result_rect.assign_subblock(hrect_ao, rf_start1, rf_end1, cf_start2, cf_end2,
783  rf_start12, cf_start12);
784  }
785  }
786 
787  result_ = ResultFactory::create(brabasis->so_matrixkit(),braaodim,ketaodim);
788  ResultFactory::convert(result_,result_rect);
789 
790  }
791 
792  }
793 
794  const ResultType& result() const { return result_; }
795 
796  private:
797 
798  ResultType result_;
799 
800  }; // class OneBodyFockMatrixBuilder
801 
802 } // end of namespace sc
803 
804 #endif // end of header guard
805 
806 // Local Variables:
807 // mode: c++
808 // c-file-style: "CLJ-CONDENSED"
809 // End:
sc::overlap
RefSCMatrix overlap(const OrbitalSpace &space2, const OrbitalSpace &space1, const Ref< SCMatrixKit > &kit=SCMatrixKit::default_matrixkit())
overlap(s2,s1) returns the overlap matrix between s2 and s1.
sc::TwoBodyFockMatrixBuilder::TwoBodyFockMatrixBuilder
TwoBodyFockMatrixBuilder(bool compute_F, bool compute_J, bool compute_K, const Ref< GaussianBasisSet > &brabasis, const Ref< GaussianBasisSet > &ketbasis, const Ref< GaussianBasisSet > &densitybasis, const RefSymmSCMatrix &density, const RefSymmSCMatrix &openshelldensity, const Ref< Integral > &integral, const Ref< MessageGrp > &msg, const Ref< ThreadGrp > &thr, double accuracy=1e-20)
Definition: fockbuilder.h:327
sc::RefSymmSCMatrix
The RefSymmSCMatrix class is a smart pointer to an SCSymmSCMatrix specialization.
Definition: matrix.h:265
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition: matrix.h:135
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::symmetrize
void symmetrize(const Ref< GPetiteList2 > &plist2, const Ref< Integral > &integral, const RefSymmSCMatrix &skel, const RefSymmSCMatrix &sym)
Uses plist2 to convert the "skeleton" matrix into the full matrix. Only applicable when the two basis...
sc::create
DescribedClass * create()
This is used to pass a function that make void constructor calls to the ClassDesc constructor.
Definition: class.h:102
sc::GaussianBasisSetMap::fblock_size
int fblock_size(int b) const
the number of basis functions in fblock b
sc::RefSymmSCMatrix::clone
RefSymmSCMatrix clone() const
These call the SCMatrix members of the same name after checking for references to 0.
sc::FeatureNotImplemented
This is thrown when an attempt is made to use a feature that is not yet implemented.
Definition: scexception.h:120
sc::SymmOneBodyIntIter
Iterator over symmetry unique shell pairs.
Definition: symmint.h:41
sc::convert
std::vector< double > convert(const RefDiagSCMatrix &A)
Converts RefDiagSCMatrix to std::vector<double>
sc::RefSCDimension
The RefSCDimension class is a smart pointer to an SCDimension specialization.
Definition: dim.h:152
sc::TwoBodyFockMatrixTransformBuilder
Builds the two-body part of the Fock matrix in MO basis using AO->MO transforms.
Definition: fockbuilder.h:495
sc::GaussianBasisSetMap::map_function
int map_function(int f) const
maps function f in Source to its location in Target
sc::TwoBodyFockMatrixBuilder
Builds the two-body part of the Fock matrix in AO basis using integral-direct algorithm.
Definition: fockbuilder.h:320
sc::GaussianBasisSetMap::fblock_to_function
int fblock_to_function(int b) const
returns the Source index of the first basis function in fblock b
sc::OneBodyFockMatrixBuilder
Builds the one-body part of the Fock matrix in AO basis.
Definition: fockbuilder.h:685
sc::FockDistribution
FockDistribution is a factory for constructing the desired FockDist specialization.
Definition: fockdist.h:187
sc::RefSymmSCMatrix::accumulate_transform
void accumulate_transform(const RefSCMatrix &a, const RefSymmSCMatrix &b, SCMatrix::Transform=SCMatrix::NormalTransform) const
Add a * b * a.t() to this.
sc::TwoBodyFockMatrixDFBuilder
Builds the two-body part of the Fock matrix in AO basis using DF-based algorithm.
Definition: fockbuilder.h:561
sc::Integral
The Integral abstract class acts as a factory to provide objects that compute one and two electron in...
Definition: integral.h:111
sc::RefSCMatrix::t
RefSCMatrix t() const
Return the transpose of this.
sc::GaussianBasisSetMap
A heavy-duty map from one GaussianBasisSet to another GaussianBasisSet.
Definition: gaussbas.h:640
sc::FockBuild
The FockBuild class works with the FockBuildThread class to generate Fock matrices for both closed sh...
Definition: fockbuild.h:805
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition: scexception.h:92
sc::SCMatrix::Transform
Transform
types of matrix transforms. Only real-valued matrices are assumed.
Definition: abstract.h:205
sc::HSOSHFContribution
Computes components of the Fock matrix necessary for high-spin open-shell calculations (e....
Definition: hsoshfcontrib.h:24
sc::GaussianBasisSetMap::nfblock
int nfblock() const
it is useful to organize sets of basis functions that are adjacent in Source and in Target.
sc::CLHFContribution
Computes components of the Fock matrix necessary for closed-shell calculations (i....
Definition: clhfcontrib.h:17
sc::OneBodyIntOp
Definition: obint.h:300
sc::transform
RefSCMatrix transform(const OrbitalSpace &space2, const OrbitalSpace &space1, const Ref< SCMatrixKit > &kit=SCMatrixKit::default_matrixkit())
transform(s2,s1) returns matrix that transforms s1 to s2.
sc::detail::FockMatrixType
Definition: fockbuilder.h:46
sc::RefCount
The base class for all reference counted objects.
Definition: ref.h:192
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
Generated at Sun Jan 26 2020 23:23:58 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.