MPQC  3.0.0-alpha
ccr12_info.h
1 //
2 // ccr12_info.h -- common utilities for all CC/CC-R12 methods
3 //
4 // Copyright (C) 2009 Toru Shiozaki
5 //
6 // Author: Toru Shiozaki <shiozaki.toru@gmail.com>
7 // Maintainer: TS & EFV
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 _chemistry_qc_ccr12_ccr12_info_h
29 #define _chemistry_qc_ccr12_ccr12_info_h
30 
31 #include <string>
32 #include <vector>
33 #include <util/misc/compute.h>
34 #include <util/group/memory.h>
35 #include <util/group/message.h>
36 #include <util/group/thread.h>
37 #include <chemistry/qc/basis/obint.h>
38 #include <chemistry/qc/basis/tbint.h>
39 #include <chemistry/qc/scf/scf.h>
40 #include <chemistry/qc/mbptr12/mbptr12.h>
41 #include <chemistry/qc/mbptr12/r12wfnworld.h>
42 #include <chemistry/qc/mbptr12/r12int_eval.h>
43 #include <chemistry/qc/ccr12/tensor.h>
44 #include <chemistry/qc/ccr12/tensorextrap.h>
45 
46 namespace sc {
47 
50 class CCR12_Info : virtual public RefCount {
51  protected:
52  const Ref<SCF> ref_;
53  const Ref<MemoryGrp>& mem_;
54 
55  const std::string theory_;
56  const std::string perturbative_;
57 
58  // common tensors
59  Ref<Tensor> d_f1;
60  Ref<Tensor> d_v2; bool need_w2_; bool need_w1_;
61  Ref<Tensor> d_t1; bool need_t1_;
62  Ref<Tensor> d_t2; bool need_t2_;
63  Ref<Tensor> d_gt2; bool need_gt2_;
64  Ref<Tensor> d_t3; bool need_t3_;
65  Ref<Tensor> d_t4; bool need_t4_;
66  Ref<Tensor> d_vr2; bool need_VpA_;// V^gg_ii ; sometimes V with CABS indices is not needed.
67  Ref<Tensor> d_vd2; // V^ii_gg
68  Ref<Tensor> d_bs2; // B^ii_ii
69  Ref<Tensor> d_xs2; // X^ii_ii
70  Ref<Tensor> d_ps2; // P^ii_ii
71  Ref<Tensor> d_fr2; bool need_FAA_; bool need_F_; // F12^AA_ii --- sometimes F with two CABS indices is not needed.
72  Ref<Tensor> d_fd2; // F12^ii_AA
73 
74  Ref<Tensor> d_vd2_gen;
75 
76  Ref<Tensor> d_qy; // F12 * gt2
77  Ref<Tensor> d_qx; // F12 * gt2 with two CABS (for CCSD-R12)
78  Ref<Tensor> d_ly; // gl2 * F12^dagger
79  Ref<Tensor> d_lx; // gl2 * F12^dagger with two CABS (for CCSD-R12)
80 
81  Ref<Tensor> d_lambda1;
82  Ref<Tensor> d_lambda2;
83  Ref<Tensor> d_glambda2;
84  Ref<Tensor> d_lambda3;
85  Ref<Tensor> d_lambda4;
86 
87  // computes source integrals
88  void compute_source_integrals_rhf();
89  void compute_source_integrals_uhf();
90  void compute_source_integrals_rhf_r12();
91 
92  // this is the OrbitalSpace representing the correlated MO space assumed by SMITH
93  Ref<OrbitalSpace> corr_space_; // full space
94  Ref<OrbitalSpace> aobs_space_; // alpha spin only, no CABS
95  Ref<OrbitalSpace> bobs_space_; // beta spin only, no CABS
96  void compute_corr_space();
97 
98  // Fgg
99  RefSCMatrix F_[NSpinCases1];
100  // FgA and FAg
101  RefSCMatrix FpA_[NSpinCases1], FAp_[NSpinCases1];
102  // source <pp|ERI|pp> integrals
103  Ref<DistArray4> pppp_acc_[NSpinCases2];
104  // this computes R12 intermediates
105  Ref<R12IntEval> r12int_eval_;
106  // V_iigg
107  RefSCMatrix Vgg_[NSpinCases2];
108  // source <pp|ERI|pA> integrals
109  Ref<DistArray4> pppA_acc_[NSpinCases2];
110  // source <ii|F12|aA> integrals
111  Ref<DistArray4> iiaA_acc_[NSpinCases2];
112 
113 
114  long irrep_f_;
115  long irrep_v_;
116  long irrep_t_;
117  long irrep_e_;
118  long irrep_y_;
119 
121  void offset_f1();
122  void offset_v2();
123  void offset_vr2();
124  void offset_vd2();
125  void offset_fr2();
126  void offset_fd2();
127  void offset_vd2_gen(bool need_cabs, bool need_xx);
128 
129  // tiling data
130  int maxtilesize_;
131  std::vector<long> sym_;
132  std::vector<long> range_;
133  std::vector<long> spin_;
134  std::vector<long> alpha_;
135  std::vector<long> offset_;
136  long noab_, nvab_, ncab_;
137 
138  // global data
139  bool restricted_;
140  int nirrep_;
141  bool fixed_;
142 
143  // orbital data
144  int nfzc_, nfzv_;
145  int naoa_, naob_;
146  int nava_, navb_;
147  int nria_, nrib_;
148  std::vector<unsigned int> mosyma_,mosymb_;
149  std::vector<unsigned int> cabssyma_,cabssymb_;
150  const Ref<R12WavefunctionWorld>& r12world_;
151  size_t workmemsize_;
152 
153  // orbital data after sorted
154  std::vector<double> orbital_evl_sorted_;
155  std::vector<long> orbital_spin_sorted_;
156  std::vector<long> orbital_sym_sorted_;
157  std::vector<long> momap_;
158 
159  // functions for initialization
160  void set_naocc(int i, int j){naoa_=i; naob_=j;};
161  void set_navir(int i, int j){nava_=i; navb_=j;};
162  void set_ncabs(int i, int j){nria_=i; nrib_=j;};
163  void set_mosym(std::vector<unsigned int> i, std::vector<unsigned int> j){mosyma_=i; mosymb_=j;};
164  void set_cabssym(std::vector<unsigned int> i, std::vector<unsigned int> j){cabssyma_=i; cabssymb_=j;};
165  void set_fixed(bool fixed){fixed_=fixed;};
166  void determine_tilesizes();
167  int determine_tilesize_each(int, int, int, std::vector<unsigned int>&, bool, int&, int);
168  void determine_maxtilesize(double);
169  void print_tile_info();
170  void needs();
171  void orbital_energies();
172 
174  void transpose( double* dest,const double* source,const int dim1,const int dim2,const double factor);
175  void transpose_1(double* dest,const double* source,const int dim1,const int dim2);
176 
178  void form_ad(Ref<Tensor>& out);
179  void form_ca(const Ref<Tensor>&, const Ref<Tensor>&, Ref<Tensor>&);
180  void form_adt(const Ref<Tensor>&, const Ref<Tensor>&, Ref<Tensor>&);
181 
184  RefSymmSCMatrix X_;
185  void retrieve_B_and_X_ii();
186  RefSymmSCMatrix B_ip_;
187  RefSymmSCMatrix X_ip_;
188  void retrieve_B_and_X_ip();
189 
190  RefDiagSCMatrix bdiag_;
191  RefSCMatrix lmatrix_;
192 
193  public:
195  const Ref<SCF>,int,int,int,long,long,int,int,
196  std::string,std::string,int);
197  ~CCR12_Info();
198 
199  void print(std::ostream&);
200 
201  double magnetic_moment() const { return ref_->magnetic_moment(); }
202 
204  const Ref<R12IntEval>& r12eval() const { return r12int_eval_; }
205  const Ref<R12WavefunctionWorld>& r12world() { return r12world_; }
206  const Ref<SCF> ref(){return ref_;};
207  bool fixed() const {return fixed_;};
208  long nirrep()const {return (long)nirrep_;};
209  int nirrep_int() const {return nirrep_;};
210  int nfzc() const {return nfzc_;};
211  int nfzv() const {return nfzv_;};
212  int naoa() const {return naoa_;};
213  int nava() const {return nava_;};
214  int naob() const {return naob_;};
215  int navb() const {return navb_;};
216  int nria() const {return nria_;};
217  int nrib() const {return nrib_;};
218  int mosyma(int i) const {return static_cast<int>(mosyma_[i]);};
219  int mosymb(int i) const {return static_cast<int>(mosymb_[i]);};
220  int cabssyma(int i) const {return static_cast<int>(cabssyma_[i]);};
221  int cabssymb(int i) const {return static_cast<int>(cabssymb_[i]);};
222 
223 
225  long noab() const {return noab_;};
226  long nvab() const {return nvab_;};
227  long ncab() const {return ncab_;};
228  long nab() const {return noab_+nvab_+ncab_;};
229  bool restricted() const {return restricted_;};
230 
231  long maxtilesize() const { return (long)maxtilesize_; };
232  long get_sym(long tile) const {return sym_[tile];};
233  long get_range(long tile) const {return range_[tile];};
234  long get_spin(long tile) const {return spin_[tile];};
235  long get_offset(long tile)const {return offset_[tile];};
236  long get_alpha(long tile) const {return alpha_[tile];};
237  double get_orb_energy(long orb) const {return orbital_evl_sorted_[orb];};
238 
239  Ref<Tensor> f1() const {return d_f1;};
240  Ref<Tensor> v2() const {return d_v2;};
241  Ref<Tensor> t1() const {return d_t1;};
242  Ref<Tensor> t2() const {return d_t2;};
243  Ref<Tensor> gt2()const {return d_gt2;};
244  Ref<Tensor> t3() const {return d_t3;};
245  Ref<Tensor> t4() const {return d_t4;};
246 
247  Ref<Tensor> fr2() const {return d_fr2;};
248  Ref<Tensor> fd2() const {return d_fd2;};
249  Ref<Tensor> vr2() const {return d_vr2;};
250  Ref<Tensor> vd2() const {return d_vd2;};
251  Ref<Tensor> xs2() const {return d_xs2;};
252  Ref<Tensor> bs2() const {return d_bs2;};
253  Ref<Tensor> ps2() const {return d_ps2;};
254 
255  Ref<Tensor> vd2_gen() const {return d_vd2_gen; };
256 
257  Ref<Tensor> qy() const {return d_qy;};
258  Ref<Tensor> qx() const {return d_qx;};
259  Ref<Tensor> ly() const {return d_ly;};
260  Ref<Tensor> lx() const {return d_lx;};
261 
262  Ref<Tensor> lambda1() const {return d_lambda1;};
263  Ref<Tensor> lambda2() const {return d_lambda2;};
264  Ref<Tensor> glambda2() const {return d_glambda2;};
265  Ref<Tensor> lambda3() const {return d_lambda3;};
266 
267  const Ref<MemoryGrp>& mem() const {return mem_;};
268  long irrep_f() const {return irrep_f_;};
269  long irrep_v() const {return irrep_v_;};
270  long irrep_t() const {return irrep_t_;};
271  long irrep_e() const {return irrep_e_;};
272  long irrep_y() const {return irrep_y_;};
273 
274 
276  void restricted_2(const long, const long, long&, long& ) const;
277  void restricted_4(const long, const long, const long, const long, long&, long&, long&, long&) const;
278  void restricted_6(const long, const long, const long, const long, const long, const long,
279  long&, long&, long&, long&, long&, long&) const;
280  void restricted_8(const long, const long, const long, const long, const long, const long, const long, const long,
281  long&, long&, long&, long&, long&, long&, long&, long&) const;
282  void sort_indices0(const double* a,double* b,const double facter) const {*b=facter*(*a);};
283  void sort_indices2(const double*, double*, const long, const long,
284  const int, const int, const double) const;
285  void sort_indices4(const double*, double*, const long, const long, const long, const long,
286  const int, const int, const int, const int, const double) const;
287  void sort_indices6(const double*, double*, const long, const long, const long, const long, const long, const long,
288  const int, const int, const int, const int, const int, const int, const double) const;
289  void sort_indices8(const double*, double*, const long, const long, const long, const long, const long, const long, const long, const long,
290  const int, const int, const int, const int, const int, const int, const int, const int, const double) const;
291  void sort_indices_acc6(const double*, double*, const long, const long, const long, const long, const long, const long,
292  const int, const int, const int, const int, const int, const int, const double) const;
293  void sort_indices_acc8(const double*, double*, const long, const long, const long, const long, const long, const long, const long, const long,
294  const int, const int, const int, const int, const int, const int, const int, const int, const double) const;
295  void smith_dgemm(const long,const long,const long,const double,const double*,const long,
296  const double*,const long,const double,double*,const long) const;
297 
298 
300  bool need_t1() const {return need_t1_;};
301  bool need_t2() const {return need_t2_;};
302  bool need_t3() const {return need_t3_;};
303  bool need_t4() const {return need_t4_;};
304  bool need_w1() const {return need_w1_;};
305  bool need_w2() const {return need_w2_;};
306  bool need_gt2()const {return need_gt2_;};
307  bool need_FAA()const {return need_FAA_;};
308  bool need_VpA()const {return need_VpA_;};
309 
310 
312  double get_e(const Ref<Tensor>&);
313  void jacobi_t1(const Ref<Tensor>& r1){jacobi_t1_(r1,d_t1);};
314  void jacobi_t1_(const Ref<Tensor>&,Ref<Tensor>&);
315  void jacobi_t2(const Ref<Tensor>& r2){jacobi_t2_(r2,d_t2);};
316  void jacobi_t2_(const Ref<Tensor>&,Ref<Tensor>&);
317  void jacobi_t3(const Ref<Tensor>& r3){jacobi_t3_(r3,d_t3);};
318  void jacobi_t3_(const Ref<Tensor>&,Ref<Tensor>&);
319  void jacobi_t4(const Ref<Tensor>& r4){jacobi_t4_(r4,d_t4);};
320  void jacobi_t4_(const Ref<Tensor>&,Ref<Tensor>&);
321  void jacobi_t2_and_gt2(const Ref<Tensor>& r2,const Ref<Tensor>& gr2) {jacobi_t2_and_gt2_(r2,d_t2,gr2,d_gt2);};
322  void jacobi_t2_and_gt2_(const Ref<Tensor>&,Ref<Tensor>&,const Ref<Tensor>&,Ref<Tensor>&);
323 
324  void jacobi_lambda1(const Ref<Tensor>& lr1) {jacobi_l1_(lr1,d_lambda1);};
325  void jacobi_l1_(const Ref<Tensor>&,Ref<Tensor>&);
326  void jacobi_lambda2(const Ref<Tensor>& lr2) {jacobi_l2_(lr2,d_lambda2);};
327  void jacobi_l2_(const Ref<Tensor>&,Ref<Tensor>&);
328 /* void jacobi_gl2(Ref<Tensor>&,Ref<Tensor>&);
329  void jacobi_l3(Ref<Tensor>&,Ref<Tensor>&);
330  void jacobi_l4(Ref<Tensor>&,Ref<Tensor>&); */
331 
333  double energy_lagrangian_r2(const Ref<Tensor>&) const;
334  double energy_lagrangian_r3(const Ref<Tensor>&) const;
335 
336 
338  Ref<SCExtrapData> edata(Ref<Tensor> t){return new TensorExtrapData(t);};
339  Ref<SCExtrapError> eerr(Ref<Tensor> t) {return new TensorExtrapError(t);};
340 
342  void offset_t1(Ref<Tensor>&, bool);
343  void offset_t2(Ref<Tensor>&, bool);
344  void offset_gt2(Ref<Tensor>&, bool);
345  void offset_t3(Ref<Tensor>&, bool);
346  void offset_t4(Ref<Tensor>&, bool);
347  void offset_e(Ref<Tensor>&);
348 
349  void offset_x_gen(Ref<Tensor>& t, const bool need_xx, const bool lprint = false);
350 
351  void offset_qy();
352  void update_qy();
353  void offset_qx();
354  void update_qx();
355  void offset_ly();
356  void update_ly();
357  void offset_lx();
358  void update_lx();
359 
360  void offset_l1(Ref<Tensor>&);
361  void offset_l2(Ref<Tensor>&);
362 
364  void guess_t2(Ref<Tensor>& d_t2_);
365  void guess_t2_r12(Ref<Tensor>& d_t2_, Ref<Tensor>& d_gt2_);
366  void guess_lambda1(Ref<Tensor>& d_lambda1_);
367  void guess_lambda1() {guess_lambda1(d_lambda1);};
368  void guess_lambda2(Ref<Tensor>& d_lambda2_);
369  void guess_lambda2() {guess_lambda2(d_lambda2);};
370  void guess_glambda2(Ref<Tensor>& d_glambda2_);
371  void guess_glambda2() {guess_glambda2(d_glambda2);};
372 
374  void fill_in_f1();
375  void fill_in_v2();
376  void fill_in_iiii();
377  void fill_in_vr_and_vd();
378  void fill_in_fr_and_fd();
379  void fill_in_vd2_gen(bool need_cabs, bool need_xx);
380 
381 
383  long momap(long i) const {return momap_[i];};
384 
386  void prod_iiii(const Ref<Tensor>&, const Ref<Tensor>&, Ref<Tensor>&, const bool transpose = false);
387 
389  RefSymmSCMatrix B() { return B_; };
390  RefSymmSCMatrix X() { return X_; };
391  RefSymmSCMatrix B_ip() { return B_ip_; };
392  RefSymmSCMatrix X_ip() { return X_ip_; };
393 
394  // returns shared pointers of OrbitalSpace objects
395  Ref<OrbitalSpace> corr_space() { return corr_space_; }; // full space
396 
397  // used in MP2-R12 updates etc.
398  void denom_contraction(const Ref<Tensor>&, Ref<Tensor>&);
399  void prediagon(RefDiagSCMatrix& eigvals, RefSCMatrix& eigvecs);
400  RefSCMatrix lmatrix() { return lmatrix_; };
401  RefDiagSCMatrix bdiag() { return bdiag_; };
402 };
403 
404 }
405 
406 #endif
sc::CCR12_Info::need_t1
bool need_t1() const
Constatnts used in specific (i.e. derived) CC-R12 object.
Definition: ccr12_info.h:300
sc::RefSymmSCMatrix
The RefSymmSCMatrix class is a smart pointer to an SCSymmSCMatrix specialization.
Definition: matrix.h:265
sc::CCR12_Info::energy_lagrangian_r2
double energy_lagrangian_r2(const Ref< Tensor > &) const
Functions for ddot of the tensors in an accurate way (used for energy evaluation).
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::CCR12_Info::get_e
double get_e(const Ref< Tensor > &)
Functions used in specific (i.e. derived) CC-R12 object.
sc::X
Definition: reftestx.h:32
sc::CCR12_Info::B
RefSymmSCMatrix B()
returns B and X intermediate for perturbative methods etc.
Definition: ccr12_info.h:389
sc::CCR12_Info::restricted_2
void restricted_2(const long, const long, long &, long &) const
Functions used in amplitude evaluators.
sc::CCR12_Info::transpose
void transpose(double *dest, const double *source, const int dim1, const int dim2, const double factor)
utilities for efficient sort_indices(n) to be implemented
sc::RefDiagSCMatrix
The RefDiagSCMatrix class is a smart pointer to an DiagSCMatrix specialization.
Definition: matrix.h:389
sc::CCR12_Info::momap
long momap(long i) const
utilities for fill-in routines
Definition: ccr12_info.h:383
sc::CCR12_Info::offset_f1
void offset_f1()
offset routines for const tensors
sc::CCR12_Info::B_
RefSymmSCMatrix B_
B and X intermediate in RefSymmSCMatrix; to be used in a certain class of methods.
Definition: ccr12_info.h:183
sc::CCR12_Info::r12eval
const Ref< R12IntEval > & r12eval() const
constants used in initialization
Definition: ccr12_info.h:204
sc::CCR12_Info::prod_iiii
void prod_iiii(const Ref< Tensor > &, const Ref< Tensor > &, Ref< Tensor > &, const bool transpose=false)
utilities for Lambda contribution in fixed-amp approaches
sc::CCR12_Info::offset_t1
void offset_t1(Ref< Tensor > &, bool)
offset routines
sc::CCR12_Info::guess_t2
void guess_t2(Ref< Tensor > &d_t2_)
guess routine
sc::CCR12_Info::form_ad
void form_ad(Ref< Tensor > &out)
local functions for jacobi_t2_and_gt2_ (i.e. a non-iterative MP2-R12 solver)
sc::CCR12_Info::noab
long noab() const
Constants used in amplitude evaluators.
Definition: ccr12_info.h:225
sc::RefCount
The base class for all reference counted objects.
Definition: ref.h:192
sc::CCR12_Info::fill_in_f1
void fill_in_f1()
fill-in routines
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::CCR12_Info
CCR12_Info is the compilation of members that are used in CC and CC-R12 methods.
Definition: ccr12_info.h:50

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