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