MADNESS  0.10.1
mp3.h
Go to the documentation of this file.
1 //
2 // Created by Florian Bischoff on 2/15/24.
3 //
4 
5 #ifndef MP3_H
6 #define MP3_H
7 
8 
9 #include <madness/mra/mra.h>
15 #include <algorithm>
16 #include <iomanip>
17 #include <iostream>
18 #include <madness/mra/macrotaskq.h>
19 
20 namespace madness {
21 class MP3 : public CCPotentials {
22 public:
23 
24  MP3(World& world, const std::shared_ptr<Nemo> nemo, const CCParameters& param)
25  : CCPotentials(world,nemo,param) {}
26  MP3(const CCPotentials& ops) : CCPotentials(ops) {}
27 
28  double mp3_energy_contribution(const Pairs<CCPair>& mp2pairs) const;
29 
30  /// compute the MP3 energy contribution, macrotask version
31  double mp3_energy_contribution_macrotask_driver(const Pairs<CCPair>& mp2pairs) const;
32 
33 private:
34  /// helper class for calculating the MP3 energy contributions
36 
38  public:
39  Partitioner(const std::string shape) {
42  if (shape=="triangular") dimension=1;
43  else if (shape=="square") dimension=2;
44  else {
45  std::string msg = "Unknown partitioning shape: " + shape;
46  MADNESS_EXCEPTION(msg.c_str(), 1);
47  }
48  };
49  };
50 
51  public:
52  MacroTaskMP3(const std::string shape) {
53  partitioner.reset(new Partitioner(shape));
54  }
55 
56  typedef std::tuple<
57  const std::string&,
58  const std::vector<int>&, // dummy vector of size npair or nocc for scheduling
59  const std::vector<int>&, // dummy vector of size npair or nocc for scheduling
60  const std::vector<std::vector<CCPairFunction<double,6>>>& , // all pairs ij
61  const std::vector<Function<double,3>>&,
62  const std::vector<Function<double,3>>&,
63  const CCParameters&,
64  const Molecule&,
65  const Function<double,3>&,
66  const std::vector<std::string>& > argtupleT;
67 
68  using resultT =std::shared_ptr<ScalarResult<double>>;
69 
70  resultT allocator(World& world, const argtupleT& argtuple) const {
71  return std::shared_ptr<ScalarResult<double>>(new ScalarResult<double>(world));
72  }
73 
74  resultT operator() (const std::string& diagram, // which diagram to calculate
75  const std::vector<int>& ij_vec, // dummy vector of size npair or nocc
76  const std::vector<int>& j_vec, // dummy vector of size 0 or nocc
77  const std::vector<std::vector<CCPairFunction<double,6>>>& pair_square, // all pairs ij
78  const std::vector<Function<double,3>>& mo_ket, // the orbitals
79  const std::vector<Function<double,3>>& mo_bra, // the orbitals*R2
80  const CCParameters& parameters,
81  const Molecule& molecule,
82  const Function<double,3>& Rsquare,
83  const std::vector<std::string>& argument) const {
84 
85  // the partitioner will break the input vector of pairs into single pairs
86  MADNESS_CHECK(ij_vec.size()==1);
88 
89  // nact=active occupied orbitals
90  const long nact=mo_ket.size()-parameters.freeze();
91  MADNESS_CHECK(pair_square.size()==nact*nact);
92 
93  // loop over pairs i<j (triangular) or i,j (square)
94  bool is_triangular=(j_vec.size()==0);
95  if (is_triangular) MADNESS_CHECK(partitioner->dimension==1);
96 
97  // determine the orbital indices i and j for the pair
98  int i=0, j=0;
99  if (is_triangular) {
100  // the batch index is the ij composite index [0,nact*(nact+1)-1]
101  const long ij=batch.result.begin;
102  // turn composite index ij into i and j, taking care of frozen orbitals
103  PairVectorMap tri_map=PairVectorMap::triangular_map(parameters.freeze(),mo_ket.size());
104  auto ij_to_i_and_j = [&tri_map](const int ij) { return tri_map.map[ij]; };
105  auto [ii,jj]=ij_to_i_and_j(ij);
106  i=ii;
107  j=jj;
108  } else {
109  MADNESS_CHECK(partitioner->dimension==2);
110  MADNESS_CHECK(j_vec.size()==ij_vec.size());
111  i=batch.input[0].begin+parameters.freeze();
112  j=batch.input[1].begin+parameters.freeze();
113  }
114  // print("i,j,parameters.freeze()=",i,j,parameters.freeze());
115 
116  // convert vector of vectors back to Pairs
117  PairVectorMap square_map=PairVectorMap::quadratic_map(parameters.freeze(),mo_ket.size());
118  auto clusterfunctions=Pairs<std::vector<CCPairFunction<double,6>>>::vector2pairs(pair_square,square_map);
119 
120  double result=0.0;
121  World& world=Rsquare.world();
122  if (diagram=="cd")
123  result= MP3::compute_mp3_cd(world,i,j,clusterfunctions,mo_ket,mo_bra,parameters,molecule,Rsquare,argument);
124  else if (diagram=="ef")
125  result= MP3::compute_mp3_ef(world,i,j,clusterfunctions,mo_ket,mo_bra,parameters,molecule,Rsquare,argument);
126  else if (diagram=="ghij")
127  result= MP3::compute_mp3_ghij(world,i,j,clusterfunctions,mo_ket,mo_bra,parameters,molecule,Rsquare,argument);
128  else if (diagram=="klmn")
129  result= MP3::compute_mp3_klmn(world,i,j,clusterfunctions,mo_ket,mo_bra,parameters,molecule,Rsquare,argument);
130  else {
131  std::string msg = "Unknown MP3 diagram: " + diagram;
132  MADNESS_EXCEPTION(msg.c_str(), 1);
133  }
134  auto result1=std::shared_ptr<ScalarResult<double>>(new ScalarResult<double>(world));
135  *result1=result;
136  return result1;
137 
138  };
139 
140 
141  };
142 
143 
144  double compute_mp3_cd(const Pairs<CCPair>& mp2pairs) const;
145  double compute_mp3_ef(const Pairs<CCPair>& mp2pairs) const;
146  double compute_mp3_ef_with_permutational_symmetry(const Pairs<CCPair>& mp2pairs) const;
147  double compute_mp3_ef_low_scaling(const Pairs<CCPair>& mp2pairs, const Pairs<std::vector<CCPairFunction<double,6>>> clusterfunctions) const;
148  double compute_mp3_ef_as_overlap(const Pairs<CCPair>& mp2pairs, const Pairs<std::vector<CCPairFunction<double,6>>> clusterfunctions) const;
149  double compute_mp3_ghij(const Pairs<CCPair>& mp2pairs, const Pairs<std::vector<CCPairFunction<double,6>>> clusterfunctions) const;
150  double compute_mp3_ghij_fast(const Pairs<CCPair>& mp2pairs, const Pairs<std::vector<CCPairFunction<double,6>>> clusterfunctions) const;
151  double compute_mp3_klmn(const Pairs<CCPair>& mp2pairs) const;
152  double compute_mp3_klmn_fast(const Pairs<CCPair>& mp2pairs) const;
153  double mp3_test(const Pairs<CCPair>& mp2pairs, const Pairs<std::vector<CCPairFunction<double,6>>> clusterfunctions) const;
154 
155  /// compute the cd term for single pair ij
156  static double compute_mp3_cd(World& world,
157  const long i, const long j,
158  const Pairs<std::vector<CCPairFunction<double,6>>>& pair_square,
159  const std::vector<Function<double,3>>& mo_ket,
160  const std::vector<Function<double,3>>& mo_bra,
161  const CCParameters& parameters,
162  const Molecule& molecule,
163  const Function<double,3>& Rsquare,
164  const std::vector<std::string>& argument);
165 
166  /// compute the ef term for single pair ij
167  static double compute_mp3_ef(World& world,
168  const long i, const long j,
169  const Pairs<std::vector<CCPairFunction<double,6>>>& pair_square,
170  const std::vector<Function<double,3>>& mo_ket,
171  const std::vector<Function<double,3>>& mo_bra,
172  const CCParameters& parameters,
173  const Molecule& molecule,
174  const Function<double,3>& Rsquare,
175  const std::vector<std::string>& argument);
176 
177  /// compute the ghij term for single pair ij
178  ///
179  /// the term actually scales linearly with the number of occupied orbitals i, so for all i!=j return zero
180  static double compute_mp3_ghij(World& world,
181  const long i, const long j,
182  const Pairs<std::vector<CCPairFunction<double,6>>>& pair_square,
183  const std::vector<Function<double,3>>& mo_ket,
184  const std::vector<Function<double,3>>& mo_bra,
185  const CCParameters& parameters,
186  const Molecule& molecule,
187  const Function<double,3>& Rsquare,
188  const std::vector<std::string>& argument);
189 
190  /// compute the klmn term for single pair ij
191  static double compute_mp3_klmn(World& world,
192  const long i, const long j,
193  const Pairs<std::vector<CCPairFunction<double,6>>>& pair_square,
194  const std::vector<Function<double,3>>& mo_ket,
195  const std::vector<Function<double,3>>& mo_bra,
196  const CCParameters& parameters,
197  const Molecule& molecule,
198  const Function<double,3>& Rsquare,
199  const std::vector<std::string>& argument);
200 
201 };
202 }
203 
204 
205 #endif //MP3_H
long size() const
Definition: macrotaskpartitioner.h:79
long begin
Definition: macrotaskpartitioner.h:67
Batch_1D result
Definition: macrotaskpartitioner.h:136
std::vector< Batch_1D > input
Definition: macrotaskpartitioner.h:135
a 6D function, either in full or low rank form, possibly including an 2-particle function
Definition: ccpairfunction.h:373
World & world() const
Returns the world.
Definition: mra.h:648
Partitioner(const std::string shape)
Definition: mp3.h:39
helper class for calculating the MP3 energy contributions
Definition: mp3.h:35
resultT operator()(const std::string &diagram, const std::vector< int > &ij_vec, const std::vector< int > &j_vec, const std::vector< std::vector< CCPairFunction< double, 6 >>> &pair_square, const std::vector< Function< double, 3 >> &mo_ket, const std::vector< Function< double, 3 >> &mo_bra, const CCParameters &parameters, const Molecule &molecule, const Function< double, 3 > &Rsquare, const std::vector< std::string > &argument) const
Definition: mp3.h:74
std::tuple< const std::string &, const std::vector< int > &, const std::vector< int > &, const std::vector< std::vector< CCPairFunction< double, 6 > > > &, const std::vector< Function< double, 3 > > &, const std::vector< Function< double, 3 > > &, const CCParameters &, const Molecule &, const Function< double, 3 > &, const std::vector< std::string > & > argtupleT
Definition: mp3.h:66
std::shared_ptr< ScalarResult< double > > resultT
Definition: mp3.h:68
resultT allocator(World &world, const argtupleT &argtuple) const
Definition: mp3.h:70
MacroTaskMP3(const std::string shape)
Definition: mp3.h:52
Definition: mp3.h:21
double compute_mp3_cd(const Pairs< CCPair > &mp2pairs) const
Definition: mp3.cc:8
double compute_mp3_ef_with_permutational_symmetry(const Pairs< CCPair > &mp2pairs) const
compute the EF term of the MP3 energy with permutational symmetry
Definition: mp3.cc:168
double compute_mp3_klmn_fast(const Pairs< CCPair > &mp2pairs) const
Definition: mp3.cc:467
double compute_mp3_klmn(const Pairs< CCPair > &mp2pairs) const
Definition: mp3.cc:537
double compute_mp3_ghij_fast(const Pairs< CCPair > &mp2pairs, const Pairs< std::vector< CCPairFunction< double, 6 >>> clusterfunctions) const
Definition: mp3.cc:400
MP3(World &world, const std::shared_ptr< Nemo > nemo, const CCParameters &param)
Definition: mp3.h:24
double mp3_energy_contribution_macrotask_driver(const Pairs< CCPair > &mp2pairs) const
compute the MP3 energy contribution, macrotask version
Definition: mp3.cc:928
MP3(const CCPotentials &ops)
Definition: mp3.h:26
double mp3_test(const Pairs< CCPair > &mp2pairs, const Pairs< std::vector< CCPairFunction< double, 6 >>> clusterfunctions) const
Definition: mp3.cc:613
double compute_mp3_ef_low_scaling(const Pairs< CCPair > &mp2pairs, const Pairs< std::vector< CCPairFunction< double, 6 >>> clusterfunctions) const
Definition: mp3.cc:253
double compute_mp3_ghij(const Pairs< CCPair > &mp2pairs, const Pairs< std::vector< CCPairFunction< double, 6 >>> clusterfunctions) const
Definition: mp3.cc:317
double compute_mp3_ef_as_overlap(const Pairs< CCPair > &mp2pairs, const Pairs< std::vector< CCPairFunction< double, 6 >>> clusterfunctions) const
double mp3_energy_contribution(const Pairs< CCPair > &mp2pairs) const
Definition: mp3.cc:886
double compute_mp3_ef(const Pairs< CCPair > &mp2pairs) const
Definition: mp3.cc:45
Definition: macrotaskq.h:716
Batch batch
Definition: macrotaskq.h:718
std::shared_ptr< MacroTaskPartitioner > partitioner
Definition: macrotaskq.h:720
partition one (two) vectors into 1D (2D) batches.
Definition: macrotaskpartitioner.h:190
std::size_t dimension
partition one or two vectors
Definition: macrotaskpartitioner.h:199
std::size_t min_batch_size
minimum batch size
Definition: macrotaskpartitioner.h:195
std::size_t max_batch_size
maximum batch size (for memory management)
Definition: macrotaskpartitioner.h:196
Definition: molecule.h:124
helper class for returning the result of a task, which is not a madness Function, but a simple scalar
Definition: macrotaskq.h:43
A parallel world class.
Definition: world.h:132
Declares the macrotaskq and MacroTaskBase classes.
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition: madness_exception.h:190
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
Main include file for MADNESS and defines Function interface.
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
Definition: CCStructures.h:199
long freeze() const
Definition: CCStructures.h:355
Definition: CCStructures.h:381
static PairVectorMap triangular_map(const int nfreeze, const int nocc)
Definition: CCStructures.h:387
std::vector< std::pair< int, int > > map
maps pair index (i,j) to vector index k
Definition: CCStructures.h:383
static PairVectorMap quadratic_map(const int nfreeze, const int nocc)
Definition: CCStructures.h:397
Definition: CCStructures.h:421
InputParameters param
Definition: tdse.cc:203
double ij(int64_t i, int64_t j)
Definition: test_distributed_matrix.cc:8
static Molecule molecule
Definition: testperiodicdft.cc:38