MADNESS 0.10.1
CCPotentials.h
Go to the documentation of this file.
1/*
2 * CCPotentials.h
3 *
4 * Created on: 4 Jan 2017
5 * Author: kottmanj
6 */
7
8#ifndef SRC_APPS_CHEM_CCPOTENTIALS_H_
9#define SRC_APPS_CHEM_CCPOTENTIALS_H_
10
11#include<madness/chem/ccpairfunction.h>
12#include<madness/chem/CCStructures.h>
13#include<madness/chem/nemo.h>
14#include<madness/chem/projector.h>
15#include<madness/chem/SCFOperators.h>
16#include<madness/chem/electronic_correlation_factor.h>
17
18namespace madness {
19/// Class which calculates all types of CC2 Potentials
20class CCPotentials {
21public:
22 CCPotentials(World& world_, const std::shared_ptr<const Nemo> nemo, const CCParameters& param);
23
24 void reset_nemo(const std::shared_ptr<const Nemo> nemo) {
25 nemo_ = nemo;
26 mo_ket_ = (make_mo_ket(*nemo));
27 mo_bra_ = (make_mo_bra(*nemo));
28 orbital_energies_ = init_orbital_energies(*nemo);
29 };
30
31 Info update_info(const CCParameters& parameters, const std::shared_ptr<const Nemo> nemo) const {
32 Info info;
33 info.mo_bra = mo_bra().get_vecfunction();
34 info.mo_ket = mo_ket().get_vecfunction();
35 info.molecular_coordinates = nemo->get_calc()->molecule.get_all_coords_vec();
36 info.parameters = parameters;
37 info.R_square = nemo->R_square;
38 info.R = nemo->R;
39 info.U1 = nemo->ncf->U1vec();
40 info.U2 = nemo->ncf->U2();
41 info.intermediate_potentials = get_potentials;
42 info.orbital_energies = orbital_energies_;
43 info.fock=nemo->compute_fock_matrix(nemo->get_calc()->amo, nemo->get_calc()->aocc);
44 return info;
45 }
46
47 virtual
48 ~CCPotentials() {};
49
50 /// forms the regularized functions from Q and Qt Ansatz for CIS(D) where tau=0 and t=mo so that Qt=Q
51 void test_pair_consistency(const CCPairFunction<double, 6>& u, const size_t i, const size_t j,
52 const CC_vecfunction& x) const;
53
54 bool test_compare_pairs(const CCPair& pair1, const CCPair& pair2) const;
55
56 void test_pairs();
57
58 void test_singles_potential(Info& info) const;
59
60 void test();
61
62 /// Returns the full 6D function from a CCPair structure (only recommended for debug purposes)
63 real_function_6d make_6D_pair(const CCPair& pair) const;
64
65 /// Function to load a function from disc
66 /// @param[in] f the function which will be loaded
67 /// @param[in] name of the file in which the function was stored
68 /// @param do_print
69 /// @return true or false depending on if the data was found on disc
70 template <typename T, size_t NDIM>
71 bool load_function(Function<T, NDIM>& f, const std::string name, bool do_print) const {
72 bool exists = archive::ParallelInputArchive<
73 archive::BinaryFstreamInputArchive>::exists(world, name.c_str());
74 if (exists) {
75 if ((world.rank() == 0) and do_print) print("loading function", name);
76 archive::ParallelInputArchive<archive::BinaryFstreamInputArchive> ar(world, name.c_str());
77 ar & f;
78 if (do_print) f.print_size(name);
79 if (f.is_compressed()) {
80 if (world.rank()==0 and do_print) print("function is compressed -- reconstructing");
81 f.change_tree_state(reconstructed);
82 if (do_print) f.print_size(name+" reconstructed");
83 save(f, name);
84 }
85 f.set_thresh(FunctionDefaults<NDIM>::get_thresh());
86 f.truncate();
87 f.print_size(name);
88 return true;
89 } else {
90 if ((world.rank()==0) and do_print) print("could not find function",name);
91 }
92 return false;
93 }
94
95 /// Plotting (convenience)
96 void plot(const vector_real_function_3d& f, const std::string& msg) const;
97
98 /// Plotting (convenience)
99 void plot(const real_function_3d& f, const std::string& msg, const bool doprint = true) const;
100
101 /// print size of a function
102 template <size_t NDIM>
103 void print_size(const Function<double, NDIM>& f, const std::string& msg, const bool print = true) const {
104 if (print) f.print_size(msg);
105 }
106
107 /// returns a vector of all orbital energies
108 std::vector<double> get_orbital_energies() const { return orbital_energies_; }
109
110 /// returns epsilon_i + epsilon_j (needed for bsh operator of pairs)
111 double get_epsilon(const size_t i, const size_t j) const {
112 return orbital_energies_[i] + orbital_energies_[j];
113 }
114
115 /// returns epsilon_i + epsilon_j (needed for bsh operator of pairs)
116 static double get_epsilon(const size_t i, const size_t j, const Info& info) {
117 return info.orbital_energies[i] + info.orbital_energies[j];
118 }
119
120 /// returns a vector of all active mos without nuclear correlation factor (nemos)
121 vector_real_function_3d get_active_mo_ket() const {
122 vector_real_function_3d result;
123 for (size_t i = parameters.freeze(); i < mo_ket_.size(); i++) result.push_back(mo_ket_(i).function);
124 return result;
125 }
126
127 /// returns a vector of all active mos multiplied with the squared nuclear currelation factor: mo_bra = R^2*mo_ket
128 vector_real_function_3d get_active_mo_bra() const {
129 vector_real_function_3d result;
130 for (size_t i = parameters.freeze(); i < mo_bra_.size(); i++) result.push_back(mo_bra_(i).function);
131 return result;
132 }
133
134 /// get the corresponding mo bra vectors to a ket vector
135 vector_real_function_3d get_mo_bra(const CC_vecfunction& ket) const {
136 vector_real_function_3d result;
137 for (const auto& ktmp : ket.functions) {
138 result.push_back(mo_bra_(ktmp.first).function);
139 }
140 return result;
141 }
142
143 /// returns a specific mo
144 CCFunction<double, 3> mo_ket(const size_t& i) const {
145 return mo_ket_(i);
146 }
147
148 // returns constant mo_ket_
149 CC_vecfunction mo_ket() const {
150 return mo_ket_;
151 }
152
153 /// returns a specific mo multiplied with the squared nuclear correlation factor
154 CCFunction<double, 3> mo_bra(const size_t& i) const {
155 return mo_bra_(i);
156 }
157
158 // returns constant mo_bra_
159 CC_vecfunction mo_bra() const {
160 return mo_bra_;
161 }
162
163 /// make bra element: R^2|t>
164 vector_real_function_3d make_bra(const CC_vecfunction& t) const {
165 return make_bra(t.get_vecfunction());
166 }
167
168 vector_real_function_3d make_bra(const vector_real_function_3d& t) const {
169 vector_real_function_3d bra = mul(world, nemo_->ncf->square(), t);
170 truncate(world, bra);
171 return bra;
172 }
173
174 /// makes the t intermediates
175 /// t_i = mo_ket_(i) + factor*tau(i)
176 /// if factor!=1 then we can not use intermediates and set the type to UNDEFINED
177 CC_vecfunction make_t_intermediate(const CC_vecfunction& tau, const CCParameters& parameters) const;
178
179 /// makes the t intermediates
180 /// t_i = mo_ket_(i) + factor*tau(i)
181 /// if the core is frozen the core ti will just be mo_ket_
182 CC_vecfunction make_full_t_intermediate(const CC_vecfunction& tau) const;
183
184 /// makes the t intermediates
185 /// t_i = mo_ket_(i) + tau(i)
186 /// if the core is frozen the core ti will just be mo_ket_
187 static CC_vecfunction make_full_t_intermediate(const CC_vecfunction& tau, const Info& info);
188
189 /// makes the t intermediates
190
191 /// t_i = mo_ket_(i) + tau(i)
192 /// skip frozen orbitals
193 static CC_vecfunction make_active_t_intermediate(const CC_vecfunction& tau, const Info& info);
194
195 /// makes the t intermediates
196 /// t_i = mo_ket_(i) + tau
197 /// i = tau.i
198 CCFunction<double, 3> make_t_intermediate(const CCFunction<double, 3>& tau) const;
199
200private:
201 /// Helper function to initialize the const mo_bra and ket elements adn orbital energies
202 CC_vecfunction
203 make_mo_bra(const Nemo& nemo) const;
204
205 /// Helper function to initialize the const mo_bra and ket elements
206 CC_vecfunction
207 make_mo_ket(const Nemo& nemo) const;
208
209 /// Transfers Eigenvalues of HF MOs from Nemo (Tensor Format) to std::vector<double> format
210 std::vector<double>
211 init_orbital_energies(const Nemo& nemo) const;
212
213public:
214 /// return the regularized MP2 ansatz: |\tau_ij> = |u_ij> + Q12 f12 |ij>
215 static CCPair make_pair_mp2(World& world, const real_function_6d& u, const size_t i, const size_t j,
216 const Info& info, bool compute_Q12_f12_ij);
217
218 /// return the regularized CC2 ansatz: |\tau_ij> = |u_ij> + Q12t f12 |t_i t_j>
219 static CCPair make_pair_cc2(World& world, const real_function_6d& u,
220 const CC_vecfunction& gs_singles, const size_t i, const size_t j, const Info& info, const bool compute_Q12_f12_ij);
221
222 /// return the regularized CC2 ansatz: |x_ij> = |u_ij> + Q12t f12 |t_i t_j> + ?????
223 static CCPair make_pair_lrcc2(World& world, const CalcType& ctype,
224 const real_function_6d& u, const CC_vecfunction& gs_singles,
225 const CC_vecfunction& ex_singles, const size_t i, const size_t j, const Info& info,
226 const bool compute_Q12_f12);
227
228 // Pair functions
229
230 /// creates pair functions for ground states
231 /// the full pair-functions is: pair = u + Qt12f12|titj>
232 /// this is decomposed in different types of CCPairFunction
233 /// the CC_Pair2 class holds an vector of CCPairFunction
234 /// we use: Qt12 = 1 - Ot1 - Ot2 + Ot12 = 1 - Ot1(1-0.5*Ot2) - Ot2(1-0.5O*t1)
235 /// depending on parameters.QtAnsatz we use the t functions also in the projector
236 /// functions[0] = u , the pure 6D part
237 /// functions[1] = f12|titj>
238 /// functions[2] = -Ot1(1-0.5*Ot2)f12|titj>
239 /// functions[3] = -(1-0.5*Ot1)Ot2f12|titj>
240 CCPair
241 make_pair_gs(const real_function_6d& u, const CC_vecfunction& tau, const size_t i, const size_t j) const;
242
243 /// creates pair functions for Excited states
244 /// the full pair-functions is: pair = u + Qt12f12|xitj> + Qt12f12|tixj> - OxQtf12|titj> - QtOxf12|titj>
245 /// the old Ansatz (not QtAnsatz) was: pair = u + Q12f12|xitj> + Q12f12|tixj>
246 /// this is decomposed in different types of CCPairFunction
247 /// the CC_Pair2 class holds an vector of CCPairFunction
248 /// we use: Qt12 = 1 - Ot1 - Ot2 + Ot12 = 1 - Ot1(1-0.5*Ot2) - Ot2(1-0.5O*t1)
249 /// functions[0] = u , the pure 6D part
250 /// functions[1] = f12|xitj>
251 /// functions[2] = f12|tixj>
252 /// functions[3] = -Ot1(1-0.5*Ot2)f12|xitj>
253 /// functions[4] = -(1-0.5*Ot1)Ot2f12|xitj>
254 /// functions[5] = -Ot1(1-0.5*Ot2)f12|tixj>
255 /// functions[6] = -(1-0.5*Ot1)Ot2f12|tixj>
256 /// functions[7] = -OxQtf12|titj> only for QtAnsatz
257 /// functions[8] = -QtOxf12|titj> only for QtAnsatz
258 CCPair
259 make_pair_ex(const real_function_6d& u, const CC_vecfunction& tau, const CC_vecfunction& x, const size_t i,
260 const size_t j, const CalcType ctype) const;
261
262 // Energy operations
263
264 /// Compute pair correlation energies of MP2 or CC2 Ground State
265 // Off diagonal pair energies are multiplied with 2.0 to acount for their permuted partners
266 /// @param[in] u the Pair_function
267 /// @param[in] singles the Singles (for MP2 give empty function) for the energy contribution over disconnected doubles
268 /// @param[out] 2*<ij|g|u> - <ji|g|u> , where i and j are determined by u (see CC_Pair class)
269 static double
270 compute_pair_correlation_energy(World& world,
271 const CCPair& u,
272 const CC_vecfunction& singles,
273 const Info& info);
274
275 /// Compute CC2 correlation energy
276 /// @param[in] The Pair_function
277 /// @param[out] \sum_{ij} 2*<ij|g|u> - <ji|g|u> + 2*<ij|g|\tau_i\tau_j> - <ji|g|\tau_i\tau_j> , where i and j are determined by u (see CC_Pair class)
278 /// since we do not compute all pairs (symmetry reasons) the off diagonal pair energies are conted twice
279 /// the cc2 pair functions are dependent on the doubles (see CC_Pair structure, and make_pair function) so make shure they are updated
280 /// @param world
281 /// @param info
282 static double
283 compute_cc2_correlation_energy(World& world, const CC_vecfunction& singles, const Pairs<CCPair>& doubles,
284 const Info& info, const std::string msg="");
285
286
287 static double
288 compute_kinetic_energy(World& world, const vector_real_function_3d& xbra, const vector_real_function_3d& xket);
289
290 /// compute the expectation value excitation energy using the CIS/CCS/CC2 singles
291 static double
292 compute_cis_expectation_value(World& world, const CC_vecfunction& x,
293 const vector_real_function_3d& V, const bool print, const Info& info);
294
295 /// Something like a pair energy for CIS(D)/LRCC2 to estimate energy convergence
296 /// calculates the response part of s2b and s2c which are independent of the mp2 amplitudes
297 static double
298 compute_excited_pair_energy(World& world, const CCPair& d, const CC_vecfunction& x, const Info& info);
299
300 /// Compute the CIS(D) Energy Correction to CIS
301 double
302 compute_cispd_energy(const CC_vecfunction& x, const Pairs<CCPair> mp2, const Pairs<CCPair> cispd) const;
303
304 /// Compute the CC2 Excitation Energy
305 double
306 compute_cc2_excitation_energy(const CC_vecfunction& stau, const CC_vecfunction& sx, const Pairs<CCPair> dtau,
307 const Pairs<CCPair> dx) const;
308
309 // 6D operations
310
311 /// The 6D Fock residue on the cusp free pair function u_{ij}(1,2) is: (2J - Kn - Un)|u_{ij}>
312 real_function_6d
313 fock_residue_6d(const CCPair& u) const;
314
315 /// Static function for the 6D Fock residue for use in macrotask
316 static madness::real_function_6d
317 fock_residue_6d_macrotask(World& world, const CCPair& u, const CCParameters& parameters,
318 const std::vector<madness::Vector<double, 3>>& all_coords_vec,
319 const std::vector<real_function_3d>& mo_ket,
320 const std::vector<real_function_3d>& mo_bra,
321 const std::vector<real_function_3d>& U1,
322 const real_function_3d& U2);
323
324 /// Static version of make_constant_part_mp2 to be called from macrotask.
325 static madness::real_function_6d
326 make_constant_part_mp2_macrotask(World& world, const CCPair& pair, const std::vector<real_function_3d>& mo_ket,
327 const std::vector<real_function_3d>& mo_bra,
328 const CCParameters& parameters, const real_function_3d& Rsquare,
329 const std::vector<real_function_3d>& U1,
330 const std::vector<std::string> argument);
331
332 /// Compute the constant part of MP2, CC2 or LR-CC2
333 ///
334 /// depending on pair.calc_type different terms are included in the constant part.
335 /// @param[in] pair the (empty) pair function, determines the terms in the constant part, contains some bookkeeping information (bsh_eps, i, j)
336 /// @param[in] gs_singles the ground-state singles for CC2 (used for the T1-transformed SO projector), may be left empty for MP2
337 /// @param[in] ex_singles the excited-state singles for CC2 (used for the T1-transformed SO projector), may be left empty for MP2 and GS-CC2
338 /// @param[in] info the Info object, containing the some basic quantities (MOs, parameters, etc)
339 /// @return the constant part of the MP2, CC2 or LR-CC2: G(Q12(g~|titj>))
340 static madness::real_function_6d
341 make_constant_part_macrotask(World& world, const CCPair& pair,
342 const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles,
343 const Info& info);
344
345
346 /// Static function to iterate the mp2 pairs from macrotask
347 static madness::real_function_6d
348 update_pair_mp2_macrotask(World& world, const CCPair& pair, const CCParameters& parameters,
349 const std::vector<madness::Vector<double, 3>>& all_coords_vec,
350 const std::vector<real_function_3d>& mo_ket,
351 const std::vector<real_function_3d>& mo_bra,
352 const std::vector<real_function_3d>& U1,
353 const real_function_3d& U2, const real_function_6d& mp2_coupling);
354
355
356 /// iterate a pair for MP2, CC2, LRCC2 on constant singles
357 static CCPair iterate_pair_macrotask(World& world,
358 const CCPair& pair, const CC_vecfunction& gs_singles,
359 const CC_vecfunction& ex_singles,
360 const real_function_6d& coupling, const Info& info, const long maxiter);
361
362
363 /// Function evaluates the consant part of the Excited state for CIS(D) if the Q Ansatz is used
364 real_function_6d
365 make_constant_part_cispd(const CCPair& u, const CC_vecfunction& x,
366 const real_convolution_6d* Gscreen = NULL) const;
367
368 /// Function evaluates the consant part of the Excited state for CIS(D) if the Qt Ansatz is used
369 real_function_6d
370 make_constant_part_cispd_Qt(const CCPair& u, const CC_vecfunction& x,
371 const real_convolution_6d* Gscreen = NULL) const;
372
373 /// Apply the Regularization potential
374 /// \f$ V_{reg} = [ U_e - [K,f12] + f12(F12-eij) ]|titj> \f$
375 /// @param[in] ti, first function in the ket, for MP2 it is the Orbital, for CC2 the relaxed Orbital t_i=\phi_i + \tau_i
376 /// @param[in] tj, second function in the ket ...
377 /// @param[in] pointer to bsh operator (in order to screen)
378 /// @param[out] the regularization potential (unprojected), see equation above
379 real_function_6d
380 apply_Vreg(const CCFunction<double, 3>& ti, const CCFunction<double, 3>& tj,
381 const real_convolution_6d* Gscreen = NULL) const;
382
383 /// Apply the Regularization potential
384 /// \f$ V_{reg} = [ U_e - [K,f12] + f12(F12-eij) + [F,Qt] ]|titj> \f$
385 /// @param[in] ti, first function in the ket, for MP2 it is the Orbital, for CC2 the relaxed Orbital t_i=\phi_i + \tau_i
386 /// @param[in] tj, second function in the ket ...
387 /// @param[in] pointer to bsh operator (in order to screen)
388 /// @param[out] the regularization potential (unprojected), see equation above
389 std::vector<CCPairFunction<double, 6>>
390 static apply_Vreg(World& world, const CCFunction<double, 3>& ti, const CCFunction<double, 3>& tj,
391 const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles,
392 const Info& info, const std::vector<std::string>& argument, const double bsh_eps);
393
394 /// Static version of apply_Vreg to be used from a macrotask. Will eventually replace former.
395 madness::real_function_6d
396 static
397 apply_Vreg_macrotask(World& world, const std::vector<real_function_3d>& mo_ket,
398 const std::vector<real_function_3d>& mo_bra,
399 const CCParameters& parameters, const real_function_3d& Rsquare,
400 const std::vector<real_function_3d>& U1, const size_t& i, const size_t& j,
401 const FuncType& x_type, const FuncType& y_type,
402 const std::vector<std::string> argument,
403 const real_convolution_6d* Gscreen = NULL);
404
405 /// evaluates: \f$ (F(1)-ei)|ti> (x) |tj> + |ti> (x) (F(2)-ej)|tj> \f$ with the help of the singles potential
406 /// singles equation is: (F-ei)|ti> = - V(ti)
407 /// response singles equation: (F-ei-omega)|xi> = - V(xi)
408 /// response: \f$ (F12-ei-ej-omega)|xitj> = (F1 - ei - omega)|xi> (x) |tj> + |xi> (x) (F2-ej)|tj> \f$
409 /// so in both cases the result will be: |V(ti),tj> + |ti,V(tj)>
410 /// @param[in] ti, first function in the ket, for MP2 it is the Orbital, for CC2 the relaxed Orbital t_i=\phi_i + \tau_i
411 /// @param[in] tj, second function in the ket ...
412 /// @param[in] pointer to bsh operator (in order to screen)
413 real_function_6d
414 apply_reduced_F1(const CCFunction<double, 3>& ti, const CCFunction<double, 3>& tj,
415 const real_convolution_6d* Gscreen = NULL) const;
416
417 /// evaluates: \f$ (F(1)-ei)|ti> (x) |tj> + |ti> (x) (F(2)-ej)|tj> \f$ with the help of the singles potential
418 /// singles equation is: (F-ei)|ti> = - V(ti)
419 /// response singles equation: (F-ei-omega)|xi> = - V(xi)
420 /// response: \f$ (F12-ei-ej-omega)|xitj> = (F1 - ei - omega)|xi> (x) |tj> + |xi> (x) (F2-ej)|tj> \f$
421 /// so in both cases the result will be: |V(ti),tj> + |ti,V(tj)>
422 /// @param[in] ti, first function in the ket, for MP2 it is the Orbital, for CC2 the relaxed Orbital t_i=\phi_i + \tau_i
423 /// @param[in] tj, second function in the ket ...
424 /// @param[in] pointer to bsh operator (in order to screen)
425 real_function_6d
426 static apply_reduced_F(World& world, const CCFunction<double, 3>& ti, const CCFunction<double, 3>& tj,
427 const Info& info, const real_convolution_6d* Gscreen = NULL);
428
429 /// Apply Ue on a tensor product of two 3d functions: Ue(1,2) |x(1)y(2)> (will be either |ij> or |\tau_i\tau_j> or mixed forms)
430 /// The Transformed electronic regularization potential (Kutzelnigg) is R_{12}^{-1} U_e R_{12} with R_{12} = R_1*R_2
431 /// It is represented as: R_{12}^{-1} U_e R_{12} = U_e + R^-1[Ue,R]
432 /// where R^-1[Ue,R] = R^-1 [[T,f],R] (see: Regularizing the molecular potential in electronic structure calculations. II. Many-body
433 /// methods, F.A.Bischoff)
434 /// The double commutator can be evaluated as follows: R^-1[[T,f],R] = -Ue_{local}(1,2)*(Un_{local}(1) - Un_{local}(2))
435 /// @param[in] x the 3D function for particle 1
436 /// @param[in] y the 3D function for particle 2
437 /// @param[in] The BSH operator to screen: Has to be in NS form, Gscreen->modified == true
438 /// @return R^-1U_eR|x,y> the transformed electronic smoothing potential applied on |x,y> :
439 real_function_6d
440 apply_transformed_Ue(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
441 const real_convolution_6d* Gscreen = NULL) const;
442
443 /// Static version of apply_transformed_Ue for the use in a macrotask.
444 /// Will eventually replace the former.
445 real_function_6d
446 static apply_transformed_Ue_macrotask(World& world, const std::vector<real_function_3d>& mo_ket,
447 const CCParameters& parameters, const real_function_3d& Rsquare,
448 const std::vector<real_function_3d>& U1, const size_t& i, const size_t& j,
449 const FuncType& x_type, const FuncType& y_type,
450 const real_convolution_6d* Gscreen = NULL);
451
452 real_function_6d
453 static apply_Ue(World& world, const CCFunction<double, 3>& phi_i, const CCFunction<double, 3>& phi_j,
454 const Info& info, const real_convolution_6d* Gscreen);
455
456
457 static real_function_6d
458 apply_KffK(World& world, const CCFunction<double, 3>& phi_i, const CCFunction<double, 3>& phi_j,
459 const Info& info, const real_convolution_6d* Gscreen);
460 static CCPairFunction<double, 6>
461 apply_commutator_F_Qt_f12(World& world, const CCFunction<double, 3>& phi_i, const CCFunction<double, 3>& phi_j,
462 const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles,
463 const Info& info, const real_convolution_6d* Gscreen);
464
465 static CCPairFunction<double, 6>
466 apply_commutator_F_dQt_f12(World& world, const CCFunction<double, 3>& phi_i, const CCFunction<double, 3>& phi_j,
467 const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles,
468 const Info& info, const real_convolution_6d* Gscreen);
469
470 /// Apply Ue on a tensor product of two 3d functions: Ue(1,2) |x(1)y(2)> (will be either |ij> or |\tau_i\tau_j> or mixed forms)
471 /// The Transformed electronic regularization potential (Kutzelnigg) is R_{12}^{-1} U_e R_{12} with R_{12} = R_1*R_2
472 /// It is represented as: R_{12}^{-1} U_e R_{12} = U_e + R^-1[Ue,R]
473 /// where R^-1[Ue,R] = R^-1 [[T,f],R] (see: Regularizing the molecular potential in electronic structure calculations. II. Many-body
474 /// methods, F.A.Bischoff)
475 /// The double commutator can be evaluated as follows: R^-1[[T,f],R] = -Ue_{local}(1,2)*(Un_{local}(1) - Un_{local}(2))
476 /// @param[in] x the 3D function for particle 1
477 /// @param[in] y the 3D function for particle 2
478 /// @param[in] Gscreen pointer to the BSH operator in order to screen, has to be in modified NS form, Gscreen->modified==true
479 /// the f12K|xy> part will be screened with the BSH while the Kf12|xy> can not be screened with the BSH operator but maybe with the coulomb
480 /// @return R^-1U_eR|x,y> the transformed electronic smoothing potential applied on |x,y> :
481 real_function_6d
482 apply_exchange_commutator(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
483 const real_convolution_6d* Gscreen = NULL) const;
484
485 real_function_6d
486 static apply_exchange_commutator_macrotask(World& world, const std::vector<real_function_3d>& mo_ket,
487 const std::vector<real_function_3d>& mo_bra,
488 const real_function_3d& Rsquare,
489 const size_t& i, const size_t& j, const CCParameters& parameters,
490 const FuncType& x_type, const FuncType& y_type,
491 const real_convolution_6d* Gscreen = NULL);
492
493 /// This applies the exchange commutator, see apply_exchange_commutator function for information
494 real_function_6d
495 apply_exchange_commutator1(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
496 const real_convolution_6d* Gscreen = NULL) const;
497
498 /// Helper Function which performs the operation \f$ <xy|g12f12|ab> \f$
499 /// @param[in] function x, if nuclear correlation is used make sure this is the correct bra function
500 /// @param[in] function y, if nuclear correlation is used make sure this is the correct bra function
501 /// @param[in] function a,
502 /// @param[in] function b,
503 double
504 make_xy_gf_ab(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y, const CCFunction<double, 3>& a,
505 const CCFunction<double, 3>& b) const;
506
507 double make_xy_ff_ab(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
508 const CCFunction<double, 3>& a, const CCFunction<double, 3>& b) const {
509 error("xy_ff_ab not yet implemented");
510 return 0.0;
511 }
512
513 /// apply the operator gf = 1/(2\gamma)*(Coulomb - 4\pi*BSH_\gamma)
514 /// works only if f = (1-exp(-\gamma*r12))/(2\gamma)
515 static real_function_3d
516 apply_gf(World& world, const real_function_3d& f, const Info& info);
517
518 /// returns <xy|op|u>
519 /// loops over every entry in the vector and accumulates results
520 /// helper function for CIS(D) energy
521 static double
522 make_xy_op_u(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
523 const CCConvolutionOperator<double, 3>& op,
524 const std::vector<CCPairFunction<double, 6>>& u);
525
526 /// returns <xy|u> for a vector of CCPairFunction
527 /// the result is accumulated for every vercotr
528 /// helper functions for CIS(D) energy
529 static double
530 make_xy_u(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
531 const std::vector<CCPairFunction<double, 6>>& u);
532
533 /// Functions which operate with the CCPairFunction structure
534 /// @param[in] function x, if nuclear correlation is used make sure this is the correct bra function
535 /// @param[in] function y, if nuclear correlation is used make sure this is the correct bra function
536 /// @param[in] CCPairFunction u,
537 static double
538 make_xy_op_u(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
539 const CCConvolutionOperator<double, 3>& op,
540 const CCPairFunction<double, 6>& u);
541
542 /// Helper Function which returns
543 /// @return <xy|op|ab>
544 /// @param[in] function x, if nuclear correlation is used make sure this is the correct bra function
545 /// @param[in] function y, if nuclear correlation is used make sure this is the correct bra function
546 /// @param[in] function a,
547 /// @param[in] function b,
548 double
549 make_xy_op_ab(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
550 const CCConvolutionOperator<double, 3>& op, const CCFunction<double, 3>& a,
551 const CCFunction<double, 3>& b) const;
552
553 /// get the correct pair function as vector of CCPairFunction functions
554 /// @param[in] The pair functions
555 /// @param[out] The demanded pair function as vector of CCPairFunction functions (includes regularization tails)
556 static std::vector<CCPairFunction<double, 6>>
557 get_pair_function(const Pairs<CCPair>& pairs, const size_t i, const size_t j);
558
559
560 /// returns <a|g12|u>_2
561 static real_function_3d
562 apply_s2b_operation(World& world, const CCFunction<double, 3>& bra, const CCPairFunction<double, 6>& u,
563 const size_t particle, const Info& info);
564
565 /// dummy to avoid confusion and for convenience
566 real_function_6d swap_particles(const real_function_6d& f) const {
567 return madness::swap_particles<double>(f);
568 }
569
570 /// swap the particles of the CCPairFunction and return a new vector of swapped functions
571 static std::vector<CCPairFunction<double, 6>> swap_particles(const std::vector<CCPairFunction<double, 6>>& f) {
572 std::vector<CCPairFunction<double, 6>> swapped;
573 for (size_t i = 0; i < f.size(); i++) swapped.push_back(f[i].swap_particles());
574 return swapped;
575 }
576
577 /// Calculate the Overlap between two CCPairFunction <f1|f2>
578 /// @param[in] 6D function 1
579 /// @param[in] 6D function 2
580 double
581 overlap(const CCPairFunction<double, 6>& f1, const CCPairFunction<double, 6>& f2) const {
582 return inner(f1, f2, nemo_->ncf->square());
583 };
584
585 /// Computes the squared norm of the pair function <x|x>
586 /// including all regularization tails
587 /// means: if x = u + f12|mn> + |ab> --> <x|x> = <u|u> + <u|f12|mn> + <u|ab> + ....
588 double
589 overlap(const CCPair& x) const;
590
591 // Projectors
592
593 /// Apply a projector to a CC_vecfunction
594 /// result: \f$ P(tau)f[i] = \sum_k tau[k]*<k|f[i]> \f$
595 /// ket state of the projector can vary
596 /// bra state are always MOs
597 vector_real_function_3d
598 apply_projector(const CC_vecfunction& f, const CC_vecfunction& ket_) const;
599
600 /// Apply Qt projector on 6D function
601 real_function_6d
602 apply_Q12t(const real_function_6d& f, const CC_vecfunction& t) const;
603
604
605 /// Apply the Q projector to a CC_vecfunction
606 /// result: \f$ 1-c*P(tau)f[i] = 1 - c*\sum_k tau[k]*<k|f[i]> \f$
607 /// ket state of the projector can vary
608 /// bra state are always MOs
609 /// the factor c is usually 1 or 1/2
610 vector_real_function_3d
611 apply_Qt(const CC_vecfunction& f, const CC_vecfunction& ket_, const double c = 1.0) const;
612
613 /// Apply the Qt projector on a CCPairFunction
614 /// works in principle like apply_Ot
615 CCPairFunction<double, 6>
616 apply_Qt(const CCPairFunction<double, 6>& f, const CC_vecfunction& t, const size_t particle,
617 const double c = 1.0) const;
618
619 /// Apply Ot projector on decomposed or op_decomposed 6D function
620 /// The function does not work with type==pure right now (not needed)
621 /// for CCPairFunction type==decomposed_ : f = \sum_k |c_k,d_k> where c and d are stored as vectors
622 /// then the result is \sum_k |a_k,b_k> with a,b stored as vectors and
623 /// \f$ a_k = t_k \f$
624 /// \f$ b_k = Ot(d_k) or vise versa for particle==2
625 /// for CCPairFunction type == op_decomposd the function si f=op|xy> and we have for particle==1
626 /// \f$ a_k = t_k \f$
627 /// \f$ b_k = <mo_k|op|x>*y \f$
628 CCPairFunction<double, 6>
629 apply_Ot(const CCPairFunction<double, 6>& f, const CC_vecfunction& t, const size_t particle) const;
630
631 /// Apply the Greens Operator to a CCPairFunction
632 /// For CCPairFunction only type pure and type decomposed is supported
633 /// for the op_decomposed type a pure function can be constructed (not needed therefore not implemented yet)
634 real_function_6d
635 apply_G(const CCPairFunction<double, 6>& u, const real_convolution_6d& G) const;
636
637 /// Apply BSH Operator and count time
638 real_function_6d apply_G(const real_function_6d& f, const real_convolution_6d& G) const {
639 CCTimer time(world, "Apply Greens Operator");
640 const real_function_6d result = G(f);
641 time.info(true, result.norm2());
642 return result;
643 }
644
645
646 // Coupled Cluster Singles potentials which work with CCPairFunction pairs
647 // notation follows Shavit & Bartlett Many-Body Methods ...
648
649 /// Calculates the CC2 singles potential for the ground state: result = Fock_residue + V
650 /// the V part is stored in the intermediate_potentials structure
651 static vector_real_function_3d
652 get_CC2_singles_potential_gs(World& world, const CC_vecfunction& singles, const Pairs<CCPair>& doubles,
653 Info& info);
654
655 /// Calculates the CCS/CIS singles potential for the excited state: result = Fock_residue + V
656 /// the V part is stored in the intermediate_potentials structure
657 /// the expectation value is calculated and updated
658 static vector_real_function_3d
659 get_CCS_potential_ex(World& world, const CC_vecfunction& x, const bool print, Info& info);
660
661 /// Calculates the CC2 singles potential for the Excited state: result = Fock_residue + V
662 /// the V part is stored in the intermediate_potentials structure
663 static vector_real_function_3d
664 get_CC2_singles_potential_ex(World& world, const CC_vecfunction& gs_singles,
665 const Pairs<CCPair>& gs_doubles, const CC_vecfunction& ex_singles,
666 const Pairs<CCPair>& response_doubles, Info& info);
667
668 /// Calculates the CC2 singles potential for the Excited state: result = Fock_residue + V
669 /// the V part is stored in the intermediate_potentials structure
670 static vector_real_function_3d
671 get_ADC2_singles_potential(World& world, const Pairs<CCPair>& gs_doubles,
672 CC_vecfunction& ex_singles, const Pairs<CCPair>& response_doubles, Info& info);
673
674 /// The potential manager for the ground state potential
675 /// CC2 singles potential parts of the ground state
676 /// Genereal function which evaluates a CC_singles potential
677 /// @param[in] Bra vectorfunction to contract the potential
678 /// @param[in] Singles of the Ground State
679 /// @param[in] Doubles of the Ground State
680 /// @param[in] Name of the potential
681 /// @param[out] the potential (without Q application)
682 /// @param world
683 double
684 potential_energy_gs(World& world, const CC_vecfunction& bra, const CC_vecfunction& singles,
685 const Pairs<CCPair>& doubles, const PotentialType& name) const;
686
687 /// The potential manager for the ground state potential
688 /// CC2 singles potential parts of the ground state
689 /// General function which evaluates a CC_singles potential
690 /// @param[in] result_index corresponds to indices of external lines in the diagram
691 /// @param[in] Singles of the Ground State
692 /// @param[in] Doubles of the Ground State
693 /// @param[in] Name of the potential
694 /// @param[out] the potential (without Q application)
695 /// @param world
696 static std::tuple<madness::vector_real_function_3d, madness::vector_real_function_3d>
697 potential_singles_gs(World& world, const std::vector<int>& result_index, const CC_vecfunction& singles,
698 const Pairs<CCPair>& doubles,
699 const PotentialType& name, const Info& info);
700
701 /// The integra manager for the excited state potential
702 /// CC2 singles potential parts of the ground state
703 /// Genereal function which evaluates a CC_singles potential
704 /// @param[in] The bra element
705 /// @param[in] Singles of the Ground State
706 /// @param[in] Doubles of the Ground State
707 /// @param[in] Singles of the Excited State
708 /// @param[in] Doubles of the Excited State
709 /// @param[in] Name of the potential
710 /// @param[out] the potential (without Q application)
711 /// @param world
712 double
713 potential_energy_ex(World& world, const CC_vecfunction& bra, const CC_vecfunction& singles_gs,
714 const Pairs<CCPair>& doubles_gs, const CC_vecfunction& singles_ex,
715 const Pairs<CCPair>& doubles_ex, const PotentialType& name) const;
716
717 /// The potential manager for the excited state potential
718 /// CC2 singles potential parts of the ground state
719 /// Genereal function which evaluates a CC_singles potential
720 /// @param[in] result_index dummy vector for partitioning the macrotasks, must have the length of the singles vector
721 /// @param[in] Singles of the Ground State
722 /// @param[in] Doubles of the Ground State
723 /// @param[in] Singles of the Excited State
724 /// @param[in] Doubles of the Excited State
725 /// @param[in] Name of the potential
726 /// @param[out] the potential (without Q application)
727 /// @param world
728 static std::tuple<madness::vector_real_function_3d, madness::vector_real_function_3d>
729 potential_singles_ex(World& world, const std::vector<int> result_index, const CC_vecfunction& singles_gs,
730 const Pairs<CCPair>& doubles_gs, const CC_vecfunction& singles_ex,
731 const Pairs<CCPair>& doubles_ex, const PotentialType& name, const Info& info);
732
733 /// The Fock operator is partitioned into F = T + Vn + R
734 /// the fock residue R= 2J-K+Un for closed shell is computed here
735 /// J_i = \sum_k <k|r12|k> |tau_i>
736 /// K_i = \sum_k <k|r12|tau_i> |k>
737 static vector_real_function_3d
738 fock_residue_closed_shell(World& world, const CC_vecfunction& singles, const Info& info);
739
740 /// the K operator runs over ALL orbitals (also the frozen ones)
741 static real_function_3d
742 K(World& world, const CCFunction<double, 3>& f, const Info& info);
743
744 /// static version of k above for access from macrotask. will eventually replace former.
745 real_function_3d
746 static K_macrotask(World& world, const std::vector<real_function_3d>& mo_ket,
747 const std::vector<real_function_3d>& mo_bra, const real_function_3d& f,
748 const CCParameters& parameters);
749
750 /// Exchange Operator on Pair function: -(K(1)+K(2))u(1,2)
751 /// if i==j in uij then the symmetry will be exploited
752 /// !!!!Prefactor (-1) is not included here!!!!
753 real_function_6d
754 K(const real_function_6d& u, const bool symmetric) const;
755
756 /// static version of k above for access from macrotask. will eventually replace former.
757 real_function_6d
758 static K_macrotask(World& world, const std::vector<real_function_3d>& mo_ket,
759 const std::vector<real_function_3d>& mo_bra,
760 const real_function_6d& u, const bool symmetric, const CCParameters& parameters);
761
762 /// Exchange Operator on Pair function: -(K(1)+K(2))u(1,2)
763 /// K(1)u(1,2) = \sum_k <k(3)|g13|u(3,2)> |k(1)>
764 /// 1. X(3,2) = bra_k(3)*u(3,2)
765 /// 2. Y(1,2) = \int X(3,2) g13 d3
766 /// 3. result = Y(1,2)*ket_k(1)
767 /// !!!!Prefactor (-1) is not included here!!!!
768 real_function_6d
769 apply_K(const real_function_6d& u, const size_t& particle) const;
770
771 /// Static version of apply_K above for access from macrotask. Will eventually replace former.
772 real_function_6d
773 static apply_K_macrotask(World& world, const std::vector<real_function_3d>& mo_ket,
774 const std::vector<real_function_3d>& mo_bra,
775 const real_function_6d& u, const size_t& particle, const CCParameters& parameters);
776
777 /// Apply the Exchange operator on a tensor product multiplied with f12
778 /// !!! Prefactor of (-1) is not inclued in K here !!!!
779 real_function_6d
780 apply_Kf(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y) const;
781
782 /// Apply fK on a tensor product of two 3D functions
783 /// fK|xy> = fK_1|xy> + fK_2|xy>
784 /// @param[in] x, the first 3D function in |xy>, structure holds index i and type (HOLE, PARTICLE, MIXED, UNDEFINED)
785 /// @param[in] y, the second 3D function in |xy> structure holds index i and type (HOLE, PARTICLE, MIXED, UNDEFINED)
786 /// @param[in] BSH operator to screen, has to be in modified NS form, Gscreen->modified()==true;
787 real_function_6d
788 apply_fK(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
789 const real_convolution_6d* Gscreen = NULL) const;
790
791 /// Creates a 6D function with the correlation factor and two given CCFunctions
792 real_function_6d
793 make_f_xy(const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
794 const real_convolution_6d* Gscreen = NULL) const;
795
796 /// Creates a 6D function with the correlation factor and two given CCFunctions
797 real_function_6d
798 static make_f_xy(World& world, const CCFunction<double, 3>& x, const CCFunction<double, 3>& y,
799 const Info& info, const real_convolution_6d* Gscreen = NULL);
800
801 real_function_6d
802 static make_f_xy_macrotask(World& world, const real_function_3d& x_ket, const real_function_3d& y_ket,
803 const real_function_3d& x_bra, const real_function_3d& y_bra,
804 const size_t& i, const size_t& j, const CCParameters& parameters,
805 const FuncType& x_type, const FuncType& y_type,
806 const real_convolution_6d* Gscreen = NULL);
807
808 /// unprojected ccs (S3c) potential
809 /// returns 2kgtk|ti> - kgti|tk>
810 /// the ccs potential: ti = ti and tk = tauk
811 static vector_real_function_3d
812 ccs_unprojected(World& world, const CC_vecfunction& ti, const CC_vecfunction& tk, const Info& info);
813
814 /// return RMS norm and max norm of residuals
815 template <typename T, std::size_t NDIM>
816 static std::pair<double, double> residual_stats(const std::vector<Function<T, NDIM>>& residual) {
817 if (residual.size() == 0) return std::make_pair(0.0, 0.0);
818 World& world = residual.front().world();
819 auto errors = norm2s(world, residual);
820 double rnorm = 0.0, maxrnorm = 0.0;
821 for (double& e : errors) {
822 maxrnorm = std::max(maxrnorm, e);
823 rnorm += e * e;
824 }
825 rnorm = sqrt(rnorm / errors.size());
826 return std::make_pair(rnorm, maxrnorm);
827 }
828
829 static void print_convergence(const std::string name, const double rmsresidual, const double maxresidual,
830 const double energy_diff, const int iteration) {
831 const std::size_t bufsize = 255;
832 char msg[bufsize];
833 std::snprintf(msg, bufsize,
834 "convergence of %s in iteration %2d at time %8.1fs: rms/max residual, energy change %.1e %.1e %.1e",
835 name.c_str(), iteration, wall_time(), rmsresidual, maxresidual,energy_diff);
836 print(msg);
837 }
838
839 // integrals from singles potentials
840
841 /// <xi|T|ti> + <xi|Vn|ti> + 2.0*<xi,k|g|ti,k> - <xi,k|g|k,ti>
842 double
843 x_s3a(const CC_vecfunction& x, const CC_vecfunction& t) const;
844
845 /// -<xi|e_i|ti>
846 double
847 x_s3b(const CC_vecfunction& x, const CC_vecfunction& t) const;
848
849 /// 2<xik|g|itk> - <xik|g|tki>
850 double
851 x_s3c(const CC_vecfunction& x, const CC_vecfunction& t) const;
852
853 /// 2<xik|g|t1it2k> - <xik|g|t2kt1i>
854 double
855 x_s5b(const CC_vecfunction& x, const CC_vecfunction& t1, const CC_vecfunction& t2) const;
856
857 /// -(2<lk|g|it1k> - <lk|g|t1ki>)* <xi|t2l>
858 double
859 x_s5c(const CC_vecfunction& x, const CC_vecfunction& t1, const CC_vecfunction& t2) const;
860
861 /// -(2<lk|g|t3i,t1k> - <lk|g|t1k,t3i>)* <xi|t2l>
862 double
863 x_s6(const CC_vecfunction& x, const CC_vecfunction& t1, const CC_vecfunction& t2,
864 const CC_vecfunction& t3) const;
865
866 /// 2.0 <xik|g|uik>- <kxi|g|uik>
867 double
868 x_s2b(const CC_vecfunction& x, const Pairs<CCPair>& u) const;
869
870 /// -( 2.0 <xi,k|(l_g_i)(2)|ukl> - <k,xi,(l_g_i)(1)|ukl> )
871 double
872 x_s2c(const CC_vecfunction& x, const Pairs<CCPair>& u) const;
873
874 /// -( <lk|g|uik> - <kl|g|uik> )*<xi|tl>
875 double
876 x_s4a(const CC_vecfunction& x, const CC_vecfunction& t, const Pairs<CCPair>& u) const;
877
878 /// -( 2.0 <xi,k|(l_g_ti)(2)|ukl> - <k,xi,(l_g_ti)(1)|ukl> )
879 double
880 x_s4b(const CC_vecfunction& x, const CC_vecfunction& t, const Pairs<CCPair>& u) const;
881
882 /// 4.0 <xil|(k_g_tk)|uil> - 2.0 <lxi|(k_g_tk)|uil> - 2.0 <xik|(l_g_tk)|uil> + <kxi|(l_g_tk)|uil>
883 double
884 x_s4c(const CC_vecfunction& x, const CC_vecfunction& t, const Pairs<CCPair>& u) const;
885
886 /// result: \sum_k( 2<k|g|uik>_2 - <k|g|uik>_1 )
887
888 /// singles are not needed explicitly but to determine if it is response or ground state
889 /// function return a vector of size 2*singles, the first half being the s2b potential, the second half
890 /// an intermediate that will be stored by the caller in the info structure
891 ///@param world
892 ///@param external_indices
893 ///@param[in] singles:CC_vecfunction fof type response or particle (depending on this the correct intermediates will be used) the functions themselves are not needed
894 ///@param[in] doubles:Pairs of CC_Pairs (GS or Response)
895 ///@param builder
896 ///@param info
897 ///@param[out] \f$ \sum_k( 2<k|g|uik>_2 - <k|g|uik>_1 ) \f$
898 /// Q-Projector is not applied, sign is correct
899 /// if the s2b potential has already been calculated it will be loaded from the intermediate_potentials structure
900 static std::tuple<madness::vector_real_function_3d, madness::vector_real_function_3d>
901 s2b(World& world, std::vector<int> external_indices, const CC_vecfunction& singles, const Pairs<CCPair>& doubles, const
902 CCPairBuilder& builder, const Info& info);
903
904 /// result: -\sum_k( <l|kgi|ukl>_2 - <l|kgi|ukl>_1)
905
906 /// singles are not needed explicitly but to determine if it is response or ground state
907 /// function return a vector of size 2*singles, the first half being the s2c potential, the second half
908 /// an intermediate that will be stored by the caller in the info structure
909 ///@param world
910 ///@param external_index
911 ///@param[in] singles:CC_vecfunction fof type response or particle (depending on this the correct intermediates will be used) the functions themselves are not needed
912 ///@param[in] doubles:Pairs of CC_Pairs (GS or Response)
913 ///@param builder
914 ///@param info
915 ///@param[out] \f$ -\sum_k( <l|kgi|ukl>_2 - <l|kgi|ukl>_1) \f$
916 /// Q-Projector is not applied, sign is correct
917 static std::tuple<vector<Function<double, 3>>, vector<Function<double, 3>>>
918 s2c(World& world, std::vector<int> external_index, const CC_vecfunction& singles, const Pairs<CCPair>& doubles, const
919 CCPairBuilder& builder, const Info& info);
920
921 /// the S4a potential can be calcualted from the S2b potential
922 /// result is \f$ s4a_i = - <l|s2b_i>*|tau_l> \f$
923 vector_real_function_3d
924 s4a_from_s2b(const vector_real_function_3d& s2b, const CC_vecfunction& singles) const;
925
926 // result: -\sum_k( <l|kgtaui|ukl>_2 - <l|kgtaui|ukl>_1) | kgtaui = <k|g|taui>
927 ///@param world
928 ///@param external_index
929 ///@param[in] singles:CC_vecfunction fof type response or particle (depending on this the correct intermediates will be used) the functions themselves are not needed
930 ///@param[in] doubles:Pairs of CC_Pairs (GS or Response)
931 ///@param info
932 ///@param[out] \f$ -( <l|kgtaui|ukl>_2 - <l|kgtaui|ukl>_1) | kgtaui = <k|g|taui> | taui=singles_i \f$
933 /// Q-Projector is not applied, sign is correct
934 static vector_real_function_3d
935 s4b(World& world, std::vector<int> external_index, const CC_vecfunction& singles, const Pairs<CCPair>& doubles,
936 const CCPairBuilder& builder, const Info& info);
937
938
939 ///@param world
940 ///@param external_index
941 ///@param[in] singles:CC_vecfunction fof type response or particle (depending on this the correct intermediates will be used) the functions themselves are not needed
942 ///@param[in] doubles:Pairs of CC_Pairs (GS or Response)
943 ///@param info
944 ///@param[out] \f$ ( 4<l|kgtauk|uil>_2 - 2<l|kgtauk|uil>_1 - 2<k|lgtauk|uil>_2 + <k|lgtauk|uil>_1 ) \f$
945 /// Q-Projector is not applied, sign is correct
946 static vector_real_function_3d
947 s4c(World& world, std::vector<int> external_index, const CC_vecfunction& singles, const Pairs<CCPair>& doubles,
948 const CCPairBuilder& builder, const Info& info);
949
950 // update the intermediates
951 void update_intermediates(const CC_vecfunction& t) {
952 g12->update_elements(mo_bra_, t);
953 // g12.sanity();
954 f12->update_elements(mo_bra_, t);
955 // f12.sanity();
956 }
957
958 /// clear stored potentials
959 /// if a response function is given only the response potentials are cleared (the GS potentials dont change anymore)
960 void clear_potentials(const CC_vecfunction& t) const {
961 if (t.type == RESPONSE) {
962 output("Clearing Response Singles-Potentials");
963 get_potentials.clear_response();
964 }
965 else {
966 output("Clearing all stored Singles-Potentials");
967 get_potentials.clear_all();
968 }
969 }
970
971public:
972 // member variables
973 /// MPI World
974 World& world;
975 /// the nemo structure which holds the nuclear correlation factor and all other functions from the reference calculation
976 std::shared_ptr<const Nemo> nemo_;
977 /// CC_Parameters: All parameters (see class description)
978 const CCParameters& parameters;
979 /// ket functions of the reference determinant
980 CC_vecfunction mo_ket_;
981 /// bra function of the reference determinant: bra = R2*ket, with R^2 beiing the squared nuclear correlation factor
982 CC_vecfunction mo_bra_;
983 /// orbital energies
984 std::vector<double> orbital_energies_;
985 /// the coulomb operator with all intermediates
986public:
987 std::shared_ptr<CCConvolutionOperator<double, 3>> g12;
988 /// the f12 operator with all intermediates
989 std::shared_ptr<CCConvolutionOperator<double, 3>> f12;
990 /// the correlation factor, holds necessary regularized potentials
991 CorrelationFactor corrfac;
992 /// Manager for stored intermediate potentials which are s2c, s2b and the whole singles potentials without fock-residue for GS and EX state
993 mutable CCIntermediatePotentials get_potentials;
994 /// POD for basis and intermediates
995 Info info;
996
997public:
998 /// Messenger structure for formated output and to store warnings
999 CCMessenger output;
1000};
1001} /* namespace madness */
1002
1003#endif /* SRC_APPS_CHEM_CCPOTENTIALS_H_ */
plot
void plot(const char *filename, const real_function_3d &f)
Definition 3dharmonic.cc:137
SCFOperators.h
Operators for the molecular HF and DFT code.
G
Definition test_derivative.cc:24
madness::Vector
A simple, fixed dimension vector.
Definition vector.h:64
f
double(* f)(const coord_3d &)
Definition derivatives.cc:54
f1
double(* f1)(const coord_3d &)
Definition derivatives.cc:55
f2
double(* f2)(const coord_3d &)
Definition derivatives.cc:56
bufsize
const std::size_t bufsize
Definition derivatives.cc:16
truncate
std::vector< Spinor > truncate(std::vector< Spinor > arg)
Definition dirac-hatom.cc:503
electronic_correlation_factor.h
inner
std::complex< double > inner(const Fcwf &psi, const Fcwf &phi)
Definition fcwf.cc:275
K
vecfuncT K(vecfuncT &ket, vecfuncT &bra, vecfuncT &vf)
maxiter
const int maxiter
Definition gygi_soltion.cc:68
u
static double u(double r, double c)
Definition he.cc:20
save
void save(World &world, const real_function_6d &f, std::string filename)
Definition helium_exact.cc:364
bsh_eps
static double bsh_eps
Definition helium_exact.cc:62
load_function
void load_function(World &world, Function< double, NDIM > &pair, const std::string name)
Definition helium_mp2.cc:407
op
Tensor< double > op(const Tensor< double > &x)
Definition kain.cc:508
print
void print(const tensorT &t)
Definition mcpfit.cc:140
madness
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
madness::norm2s
std::vector< double > norm2s(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norms of a vector of functions.
Definition vmra.h:845
madness::mul
Function< TENSOR_RESULT_TYPE(Q, T), NDIM > mul(const Q alpha, const Function< T, NDIM > &f, bool fence=true)
Returns new function equal to alpha*f(x) with optional fence.
Definition mra.h:1795
madness::swap_particles
std::enable_if_t< NDIM%2==0, Function< T, NDIM > > swap_particles(const Function< T, NDIM > &f)
swap particles 1 and 2
Definition mra.h:2414
madness::vector_real_function_3d
std::vector< real_function_3d > vector_real_function_3d
Definition functypedefs.h:94
madness::real_function_6d
Function< double, 6 > real_function_6d
Definition functypedefs.h:83
madness::print_size
void print_size(World &world, const std::vector< Function< T, NDIM > > &v, const std::string &msg="vectorfunction")
Definition vmra.h:1818
madness::wall_time
double wall_time()
Returns the wall time in seconds relative to an arbitrary origin.
Definition timers.cc:48
madness::real_function_3d
Function< double, 3 > real_function_3d
Definition functypedefs.h:80
madness::name
std::string name(const FuncType &type, const int ex=-1)
Definition ccpairfunction.h:28
b
static const double b
Definition nonlinschro.cc:119
d
static const double d
Definition nonlinschro.cc:121
a
static const double a
Definition nonlinschro.cc:118
error
void error(const char *msg, int code)
Definition oldtest.cc:57
c
static const double c
Definition relops.cc:10
V
static double V(const coordT &r)
Definition tdse.cc:288
param
InputParameters param
Definition tdse.cc:203
e
void e()
Definition test_sig.cc:75
residual
F residual(const F &f)
Definition testcomplexfunctionsolver.cc:69
test
void test()
Definition y.cc:696
Generated on Tue Nov 4 2025 01:14:56 for MADNESS by doxygen 1.9.8