MPQC 3.0.0-alpha
Loading...
Searching...
No Matches
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
43namespace 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::CLHFContribution
Computes components of the Fock matrix necessary for closed-shell calculations (i....
Definition clhfcontrib.h:17
sc::FeatureNotImplemented
This is thrown when an attempt is made to use a feature that is not yet implemented.
Definition scexception.h:120
sc::FockBuild
The FockBuild class works with the FockBuildThread class to generate Fock matrices for both closed sh...
Definition fockbuild.h:805
sc::FockDistribution
FockDistribution is a factory for constructing the desired FockDist specialization.
Definition fockdist.h:187
sc::GaussianBasisSetMap
A heavy-duty map from one GaussianBasisSet to another GaussianBasisSet.
Definition gaussbas.h:640
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::GaussianBasisSetMap::fblock_size
int fblock_size(int b) const
the number of basis functions in fblock b
sc::GaussianBasisSetMap::nfblock
int nfblock() const
it is useful to organize sets of basis functions that are adjacent in Source and in Target.
sc::GaussianBasisSetMap::map_function
int map_function(int f) const
maps function f in Source to its location in Target
sc::HSOSHFContribution
Computes components of the Fock matrix necessary for high-spin open-shell calculations (e....
Definition hsoshfcontrib.h:24
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::OneBodyFockMatrixBuilder
Builds the one-body part of the Fock matrix in AO basis.
Definition fockbuilder.h:685
sc::OneBodyIntOp
Definition obint.h:300
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition scexception.h:92
sc::RefCount
The base class for all reference counted objects.
Definition ref.h:192
sc::RefSCDimension
The RefSCDimension class is a smart pointer to an SCDimension specialization.
Definition dim.h:152
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition matrix.h:135
sc::RefSCMatrix::t
RefSCMatrix t() const
Return the transpose of this.
sc::RefSymmSCMatrix
The RefSymmSCMatrix class is a smart pointer to an SCSymmSCMatrix specialization.
Definition matrix.h:265
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::RefSymmSCMatrix::clone
RefSymmSCMatrix clone() const
These call the SCMatrix members of the same name after checking for references to 0.
sc::Ref
A template class that maintains references counts.
Definition ref.h:361
sc::SCMatrix::Transform
Transform
types of matrix transforms. Only real-valued matrices are assumed.
Definition abstract.h:205
sc::SymmOneBodyIntIter
Iterator over symmetry unique shell pairs.
Definition symmint.h:41
sc::TwoBodyFockMatrixBuilder
Builds the two-body part of the Fock matrix in AO basis using integral-direct algorithm.
Definition fockbuilder.h:320
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::TwoBodyFockMatrixDFBuilder
Builds the two-body part of the Fock matrix in AO basis using DF-based algorithm.
Definition fockbuilder.h:561
sc::TwoBodyFockMatrixTransformBuilder
Builds the two-body part of the Fock matrix in MO basis using AO->MO transforms.
Definition fockbuilder.h:495
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
Contains all MPQC code up to version 3.
Definition mpqcin.h:14
sc::convert
std::vector< double > convert(const RefDiagSCMatrix &A)
Converts RefDiagSCMatrix to std::vector<double>
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::detail::FockMatrixType
Definition fockbuilder.h:46
Generated at Wed Sep 25 2024 02:45:29 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.12.0.