Psi4
dfhelper.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * Psi4: an open-source quantum chemistry software package
5  *
6  * Copyright (c) 2007-2018 The Psi4 Developers.
7  *
8  * The copyrights for code used from other parties are included in
9  * the corresponding files.
10  *
11  * This file is part of Psi4.
12  *
13  * Psi4 is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU Lesser General Public License as published by
15  * the Free Software Foundation, version 3.
16  *
17  * Psi4 is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU Lesser General Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public License along
23  * with Psi4; if not, write to the Free Software Foundation, Inc.,
24  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25  *
26  * @END LICENSE
27  */
28 
29 #ifndef three_index_dfhelper
30 #define three_index_dfhelper
31 
32 #include "psi4/psi4-dec.h"
33 #include <psi4/libmints/typedefs.h>
34 
35 #include <map>
36 #include <list>
37 #include <vector>
38 #include <tuple>
39 #include <string>
40 
41 namespace psi {
42 
43 class BasisSet;
44 class Options;
45 class Matrix;
46 class ERISieve;
47 class TwoBodyAOInt;
48 
50  public:
51  DFHelper(std::shared_ptr<BasisSet> primary, std::shared_ptr<BasisSet> aux);
52  ~DFHelper();
53 
62  void set_method(std::string method) { method_ = method; }
63  std::string get_method() { return method_; }
64 
70  void set_nthreads(size_t nthreads) { nthreads_ = nthreads; }
71  size_t get_nthreads() { return nthreads_; }
72 
78  void set_memory(size_t doubles) { memory_ = doubles; }
79  size_t get_memory() { return memory_; }
80 
82  size_t get_AO_size() { return big_skips_[nbf_]; }
83 
85  size_t get_core_size() {
86  AO_core();
87  return required_core_size_;
88  }
89 
91  double ao_sparsity() { return (1.0 - (double)small_skips_[nbf_] / (double)(nbf_ * nbf_)); }
92 
101  void set_AO_core(bool core) { AO_core_ = core; }
102  bool get_AO_core() { return AO_core_; }
103 
110  void set_MO_core(bool core) { MO_core_ = core; }
111  bool get_MO_core() { return MO_core_; }
112 
114  void set_schwarz_cutoff(double cutoff) { cutoff_ = cutoff; }
115  double get_schwarz_cutoff() { return cutoff_; }
116 
118  void set_metric_pow(double pow) { mpower_ = pow; }
119  double get_metric_pow() { return mpower_; }
120 
126  void hold_met(bool hold) { hold_met_ = hold; }
127  bool get_hold_met() { return hold_met_; }
128 
133  void set_fitting_condition(double condition) { condition_ = condition; }
134  bool get_fitting_condition() { return condition_; }
135 
140  void set_do_wK(bool do_wK) { do_wK_ = do_wK; }
141  size_t get_do_wK() { return do_wK_; }
142 
147  void set_omega(double omega) { omega_ = omega; }
148  size_t get_omega() { return omega_; }
149 
154  void set_print_lvl(int print_lvl) { print_lvl_ = print_lvl; }
155  int get_print_lvl() { return print_lvl_; }
156 
158  void initialize();
159 
161  void prepare_sparsity();
162 
164  void print_header();
165 
171  void add_space(std::string key, SharedMatrix space);
172 
180  void add_transformation(std::string name, std::string key1, std::string key2, std::string order = "Qpq");
181 
183  void transform();
184 
185  // => Tensor IO <=
186  // many ways to access the 3-index tensors.
187 
197  void fill_tensor(std::string name, SharedMatrix M);
198  void fill_tensor(std::string name, SharedMatrix M, std::vector<size_t> a1);
199  void fill_tensor(std::string name, SharedMatrix M, std::vector<size_t> a1, std::vector<size_t> a2);
200  void fill_tensor(std::string name, SharedMatrix M, std::vector<size_t> a1, std::vector<size_t> a2,
201  std::vector<size_t> a3);
202 
209  void fill_tensor(std::string name, double* b, std::vector<size_t> a1, std::vector<size_t> a2,
210  std::vector<size_t> a3);
211  void fill_tensor(std::string name, double* b, std::vector<size_t> a1, std::vector<size_t> a2);
212  void fill_tensor(std::string name, double* b, std::vector<size_t> a1);
213  void fill_tensor(std::string name, double* b);
214 
222  SharedMatrix get_tensor(std::string name);
223  SharedMatrix get_tensor(std::string name, std::vector<size_t> a1);
224  SharedMatrix get_tensor(std::string name, std::vector<size_t> a1, std::vector<size_t> a2);
225  SharedMatrix get_tensor(std::string name, std::vector<size_t> a1, std::vector<size_t> a2, std::vector<size_t> a3);
226 
232  void add_disk_tensor(std::string key, std::tuple<size_t, size_t, size_t> dimensions);
233 
241  void write_disk_tensor(std::string name, SharedMatrix M);
242  void write_disk_tensor(std::string name, SharedMatrix M, std::vector<size_t> a1);
243  void write_disk_tensor(std::string name, SharedMatrix M, std::vector<size_t> a1, std::vector<size_t> a2);
244  void write_disk_tensor(std::string name, SharedMatrix M, std::vector<size_t> a1, std::vector<size_t> a2,
245  std::vector<size_t> a3);
246 
254  void write_disk_tensor(std::string name, double* b, std::vector<size_t> a1, std::vector<size_t> a2,
255  std::vector<size_t> a3);
256  void write_disk_tensor(std::string name, double* b, std::vector<size_t> a1, std::vector<size_t> a2);
257  void write_disk_tensor(std::string name, double* b, std::vector<size_t> a1);
258  void write_disk_tensor(std::string name, double* b);
259 
261  void transpose(std::string name, std::tuple<size_t, size_t, size_t> order);
262 
264  void clear_spaces();
265 
267  void clear_all();
268 
270  size_t get_space_size(std::string key);
271  size_t get_tensor_size(std::string key);
272  std::tuple<size_t, size_t, size_t> get_tensor_shape(std::string key);
273  size_t get_naux() { return naux_; }
274 
276  void build_JK(std::vector<SharedMatrix> Cleft, std::vector<SharedMatrix> Cright, std::vector<SharedMatrix> D,
277  std::vector<SharedMatrix> J, std::vector<SharedMatrix> K, size_t max_nocc, bool do_J, bool do_K,
278  bool do_wK, bool lr_symmetric);
279 
280  protected:
281  // => basis sets <=
282  std::shared_ptr<BasisSet> primary_;
283  std::shared_ptr<BasisSet> aux_;
284  size_t nbf_;
285  size_t naux_;
286 
287  // => memory in doubles <=
288  size_t memory_ = 256000000;
290 
291  // => internal holders <=
292  std::string method_ = "STORE";
293  bool direct_ = false;
294  bool direct_iaQ_ = false;
296  bool AO_core_ = true;
297  bool MO_core_ = false;
298  size_t nthreads_ = 1;
299  double cutoff_ = 1e-12;
300  double condition_ = 1e-12;
301  double mpower_ = -0.5;
302  bool hold_met_ = false;
303  bool built_ = false;
304  bool transformed_ = false;
305  std::pair<size_t, size_t> info_;
306  bool ordered_ = false;
307  bool do_wK_ = false;
308  double omega_;
309  bool debug_ = false;
310  bool sparsity_prepared_ = false;
311  int print_lvl_ = 1;
312 
313  // => in-core machinery <=
314  void AO_core();
315  std::unique_ptr<double[]> Ppq_;
316  std::map<double, SharedMatrix> metrics_;
317 
318  // => AO building machinery <=
319  void prepare_AO();
320  void prepare_AO_core();
321  void compute_dense_Qpq_blocking_Q(const size_t start, const size_t stop, double* Mp,
322  std::vector<std::shared_ptr<TwoBodyAOInt>> eri);
323  void compute_sparse_pQq_blocking_Q(const size_t start, const size_t stop, double* Mp,
324  std::vector<std::shared_ptr<TwoBodyAOInt>> eri);
325  void compute_sparse_pQq_blocking_p(const size_t start, const size_t stop, double* Mp,
326  std::vector<std::shared_ptr<TwoBodyAOInt>> eri);
327  void compute_sparse_pQq_blocking_p_symm(const size_t start, const size_t stop, double* Mp,
328  std::vector<std::shared_ptr<TwoBodyAOInt>> eri);
329  void contract_metric_AO_core_symm(double* Qpq, double* metp, size_t begin, size_t end);
330  void grab_AO(const size_t start, const size_t stop, double* Mp);
331 
332  // first integral transforms
333  void first_transform_pQq(size_t bsize, size_t bcount, size_t block_size, double* Mp, double* Tp, double* Bp,
334  std::vector<std::vector<double>>& C_buffers);
335 
336  // => index vectors for screened AOs <=
337  std::vector<size_t> small_skips_;
338  std::vector<size_t> big_skips_;
339  std::vector<size_t> symm_ignored_columns_;
340  std::vector<size_t> symm_small_skips_;
341  std::vector<size_t> symm_big_skips_;
342 
343  // => shell info and blocking <=
344  size_t pshells_;
345  size_t Qshells_;
346  double Qshell_max_;
347  std::vector<size_t> pshell_aggs_;
348  std::vector<size_t> Qshell_aggs_;
349  void prepare_blocking();
350 
351  // => generalized blocking <=
352  std::pair<size_t, size_t> pshell_blocks_for_AO_build(const size_t mem, size_t symm,
353  std::vector<std::pair<size_t, size_t>>& b);
354  std::pair<size_t, size_t> Qshell_blocks_for_transform(const size_t mem, size_t wtmp, size_t wfinal,
355  std::vector<std::pair<size_t, size_t>>& b);
356  void metric_contraction_blocking(std::vector<std::pair<size_t, size_t>>& steps, size_t blocking_index,
357  size_t block_sizes, size_t total_mem, size_t memory_factor, size_t memory_bump);
358 
359  // => Schwarz Screening <=
360  std::vector<size_t> schwarz_fun_mask_;
361  std::vector<size_t> schwarz_shell_mask_;
362  std::vector<size_t> schwarz_fun_count_;
363 
364  // => Coulomb metric handling <=
365  std::vector<std::pair<double, std::string>> metric_keys_;
366  void prepare_metric();
367  void prepare_metric_core();
368  double* metric_prep_core(double pow);
369  std::string return_metfile(double pow);
370  std::string compute_metric(double pow);
371 
372  // => metric operations <=
373  void contract_metric_Qpq(std::string file, double* metp, double* Mp, double* Fp, const size_t tots);
374  void contract_metric(std::string file, double* metp, double* Mp, double* Fp, const size_t tots);
375  void contract_metric_core(std::string file);
376  void contract_metric_AO(double* Mp);
377  void contract_metric_AO_core(double* Qpq, double* metp);
378 
379  // => spaces and transformation maps <=
380  std::map<std::string, std::tuple<SharedMatrix, size_t>> spaces_;
381  std::map<std::string, std::tuple<std::string, std::string, size_t>> transf_;
382  std::map<std::string, std::unique_ptr<double[]>> transf_core_;
383 
384  // => transformation machinery <=
385  std::pair<size_t, size_t> identify_order();
386  void print_order();
387  void put_transformations_Qpq(int begin, int end, int wsize, int bsize, double* Fp, int ind, bool bleft);
388  void put_transformations_pQq(int begin, int end, int block_size, int bcount, int wsize, int bsize, double* Np,
389  double* Fp, int ind, bool bleft);
390  std::vector<std::pair<std::string, size_t>> sorted_spaces_;
391  std::vector<std::string> order_;
392  std::vector<std::string> bspace_;
393  std::vector<size_t> strides_;
394 
395  // => FILE IO maintenence <=
396  typedef struct StreamStruct {
397  StreamStruct();
398  StreamStruct(std::string filename, std::string op, bool activate = true);
399  ~StreamStruct();
400 
401  FILE* get_stream(std::string op);
402  void change_stream(std::string op);
403  void close_stream();
404 
405  FILE* fp_;
406  std::string op_;
407  bool open_ = false;
408  std::string filename_;
409 
410  } Stream;
411 
412  std::map<std::string, std::shared_ptr<Stream>> file_streams_;
413  FILE* stream_check(std::string filename, std::string op);
414 
415  // => FILE IO machinery <=
416  void put_tensor(std::string file, double* b, std::pair<size_t, size_t> a1, std::pair<size_t, size_t> a2,
417  std::pair<size_t, size_t> a3, std::string op);
418 
419  void put_tensor(std::string file, double* b, const size_t start1, const size_t stop1, const size_t start2,
420  const size_t stop2, std::string op);
421  void get_tensor_(std::string file, double* b, std::pair<size_t, size_t> a1, std::pair<size_t, size_t> a2,
422  std::pair<size_t, size_t> a3);
423  void get_tensor_(std::string file, double* b, const size_t start1, const size_t stop1, const size_t start2,
424  const size_t stop2);
425  void put_tensor_AO(std::string file, double* Mp, size_t size, size_t start, std::string op);
426  void get_tensor_AO(std::string file, double* Mp, size_t size, size_t start);
427 
428  // => internal handlers for FILE IO <=
429  std::map<std::string, std::tuple<std::string, std::string>> files_;
430  std::map<std::string, std::tuple<size_t, size_t, size_t>> sizes_;
431  std::map<std::string, std::tuple<size_t, size_t, size_t>> tsizes_;
432  std::map<std::string, std::string> AO_files_;
433  std::vector<size_t> AO_file_sizes_;
434  std::vector<std::string> AO_names_;
435  std::string start_filename(std::string start);
436  void filename_maker(std::string name, size_t a0, size_t a1, size_t a2, size_t op = 0);
437  void AO_filename_maker(size_t i);
438  void check_file_key(std::string);
439  void check_file_tuple(std::string name, std::pair<size_t, size_t> t0, std::pair<size_t, size_t> t1,
440  std::pair<size_t, size_t> t2);
441  void check_matrix_size(std::string name, SharedMatrix M, std::pair<size_t, size_t> t0, std::pair<size_t, size_t> t1,
442  std::pair<size_t, size_t> t2);
443 
444  // => transpose a tensor <=
445  void transpose_core(std::string name, std::tuple<size_t, size_t, size_t> order);
446  void transpose_disk(std::string name, std::tuple<size_t, size_t, size_t> order);
447 
448  // => JK <=
449  void compute_JK(std::vector<SharedMatrix> Cleft, std::vector<SharedMatrix> Cright, std::vector<SharedMatrix> D,
450  std::vector<SharedMatrix> J, std::vector<SharedMatrix> K, size_t max_nocc, bool do_J, bool do_K,
451  bool do_wK, bool lr_symmetric);
452  void compute_D(std::vector<SharedMatrix> D, std::vector<SharedMatrix> Cleft, std::vector<SharedMatrix> Cright);
453  void compute_J(std::vector<SharedMatrix> D, std::vector<SharedMatrix> J, double* Mp, double* T1p, double* T2p,
454  std::vector<std::vector<double>>& D_buffers, size_t bcount, size_t block_size);
455  void compute_J_symm(std::vector<SharedMatrix> D, std::vector<SharedMatrix> J, double* Mp, double* T1p, double* T2p,
456  std::vector<std::vector<double>>& D_buffers, size_t bcount, size_t block_size);
457  void compute_K(std::vector<SharedMatrix> Cleft, std::vector<SharedMatrix> Cright, std::vector<SharedMatrix> K,
458  double* Tp, double* Jtmp, double* Mp, size_t bcount, size_t block_size,
459  std::vector<std::vector<double>>& C_buffers, bool lr_symmetric);
460  std::tuple<size_t, size_t> Qshell_blocks_for_JK_build(std::vector<std::pair<size_t, size_t>>& b, size_t max_nocc,
461  bool lr_symmetric);
462 
463  // => misc <=
464  void fill(double* b, size_t count, double value);
465 
466 }; // End DF Helper class
467 } // psi4 namespace
468 #endif
std::map< std::string, std::tuple< SharedMatrix, size_t > > spaces_
Definition: dfhelper.h:380
void hold_met(bool hold)
Definition: dfhelper.h:126
std::vector< size_t > small_skips_
Definition: dfhelper.h:337
double ao_sparsity()
Returns the amount of sparsity in the AO integrals.
Definition: dfhelper.h:91
std::vector< size_t > AO_file_sizes_
Definition: dfhelper.h:433
size_t pshells_
Definition: dfhelper.h:344
void set_metric_pow(double pow)
fitting metric power (defaults to -0.5)
Definition: dfhelper.h:118
double get_metric_pow()
Definition: dfhelper.h:119
int get_print_lvl()
Definition: dfhelper.h:155
std::map< double, SharedMatrix > metrics_
Definition: dfhelper.h:316
std::vector< std::string > order_
Definition: dfhelper.h:391
void set_nthreads(size_t nthreads)
Definition: dfhelper.h:70
size_t get_do_wK()
Definition: dfhelper.h:141
void set_AO_core(bool core)
Definition: dfhelper.h:101
std::string get_method()
Definition: dfhelper.h:63
std::vector< size_t > pshell_aggs_
Definition: dfhelper.h:347
double Qshell_max_
Definition: dfhelper.h:346
std::map< std::string, std::tuple< size_t, size_t, size_t > > tsizes_
Definition: dfhelper.h:431
double get_schwarz_cutoff()
Definition: dfhelper.h:115
void set_method(std::string method)
Definition: dfhelper.h:62
Definition: pointgrp.h:104
void set_do_wK(bool do_wK)
Definition: dfhelper.h:140
bool get_hold_met()
Definition: dfhelper.h:127
std::map< std::string, std::tuple< std::string, std::string > > files_
Definition: dfhelper.h:429
FILE * fp_
Definition: dfhelper.h:405
size_t Qshells_
Definition: dfhelper.h:345
void set_fitting_condition(double condition)
Definition: dfhelper.h:133
size_t get_AO_size()
Returns the number of doubles in the screened AO integrals.
Definition: dfhelper.h:82
std::vector< std::pair< std::string, size_t > > sorted_spaces_
Definition: dfhelper.h:390
void set_MO_core(bool core)
Definition: dfhelper.h:110
std::vector< std::string > bspace_
Definition: dfhelper.h:392
std::map< std::string, std::shared_ptr< Stream > > file_streams_
Definition: dfhelper.h:412
std::vector< size_t > Qshell_aggs_
Definition: dfhelper.h:348
size_t naux_
Definition: dfhelper.h:285
std::vector< size_t > schwarz_fun_count_
Definition: dfhelper.h:362
std::string filename_
Definition: dfhelper.h:408
std::unique_ptr< double[]> Ppq_
Definition: dfhelper.h:315
std::vector< size_t > big_skips_
Definition: dfhelper.h:338
std::map< std::string, std::tuple< size_t, size_t, size_t > > sizes_
Definition: dfhelper.h:430
std::string op_
Definition: dfhelper.h:406
Definition: dfhelper.h:49
Definition: dfhelper.h:396
size_t get_nthreads()
Definition: dfhelper.h:71
bool get_MO_core()
Definition: dfhelper.h:111
std::vector< size_t > strides_
Definition: dfhelper.h:393
std::map< std::string, std::unique_ptr< double[]> > transf_core_
Definition: dfhelper.h:382
size_t required_core_size_
Definition: dfhelper.h:289
std::vector< size_t > symm_ignored_columns_
Definition: dfhelper.h:339
std::vector< std::string > AO_names_
Definition: dfhelper.h:434
std::vector< size_t > symm_small_skips_
Definition: dfhelper.h:340
std::shared_ptr< BasisSet > primary_
Definition: dfhelper.h:282
size_t get_core_size()
Returns the size of the in-core version in doubles.
Definition: dfhelper.h:85
void set_schwarz_cutoff(double cutoff)
schwarz screening cutoff (defaults to 1e-12)
Definition: dfhelper.h:114
void set_memory(size_t doubles)
Definition: dfhelper.h:78
std::map< std::string, std::tuple< std::string, std::string, size_t > > transf_
Definition: dfhelper.h:381
std::shared_ptr< BasisSet > aux_
Definition: dfhelper.h:283
void set_print_lvl(int print_lvl)
Definition: dfhelper.h:154
void set_omega(double omega)
Definition: dfhelper.h:147
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:49
#define PSI_API
Definition: pragma.h:155
double omega_
Definition: dfhelper.h:308
bool get_AO_core()
Definition: dfhelper.h:102
std::vector< size_t > schwarz_shell_mask_
Definition: dfhelper.h:361
std::pair< size_t, size_t > info_
Definition: dfhelper.h:305
size_t get_omega()
Definition: dfhelper.h:148
std::vector< size_t > symm_big_skips_
Definition: dfhelper.h:341
size_t get_memory()
Definition: dfhelper.h:79
std::vector< std::pair< double, std::string > > metric_keys_
Definition: dfhelper.h:365
std::vector< size_t > schwarz_fun_mask_
Definition: dfhelper.h:360
size_t get_naux()
Definition: dfhelper.h:273
bool space(char c)
Definition: stl_string.cc:109
size_t nbf_
Definition: dfhelper.h:284
bool get_fitting_condition()
Definition: dfhelper.h:134
std::map< std::string, std::string > AO_files_
Definition: dfhelper.h:432
bool symm_compute_
Definition: dfhelper.h:295