MADNESS 0.10.1
CC2.h
Go to the documentation of this file.
1/*
2 * CC2.h
3 *
4 * Created on: Aug 17, 2015
5 * Author: kottmanj
6 */
7
8#ifndef CC2_H_
9#define CC2_H_
10
12#include<madness/chem/SCF.h>
13#include<madness/chem/nemo.h>
15#include<madness/chem/mp3.h>
17#include <madness/mra/mra.h>
18#include <madness/mra/vmra.h>
19#include <madness/mra/lbdeux.h>
20#include <madness/misc/ran.h>
21#include<madness/chem/TDHF.h>
23
24#include "BSHApply.h"
25
26namespace madness {
27
29public:
30
31 CC2(World& world_, const commandlineparser& parser, const std::shared_ptr<Nemo> nemo_)
32 : world(world_),
33 parameters(world_,parser),
34 nemo(nemo_),
37
38 output.section("CC2 Class has been initialized with the following parameters");
39 // set the threshholds
40 // Set Protocoll
41 output("Set Protocol 3D");
42 nemo_->get_calc()->set_protocol<3>(world, parameters.thresh_3D());
43 output("Set Protocol 6D");
44 nemo_->get_calc()->set_protocol<6>(world, parameters.thresh_6D());
45
48 // Make sure that k is the same in 3d and 6d functions
50 // by default SCF sets the truncate_mode to 1
53 // set initial level to 3 (important to avoid undersampling issues)
54 // the "real initial level" is then 2 since the initial level gets substracted by one if refine is true (which is the default)
63
64 tdhf.reset(new TDHF(world,parser,nemo));
65
66 }
67
68 virtual ~CC2() {}
69
70
71 double value() {
72 return value(nemo->molecule().get_all_coords());
73 }
74
75 double value(const Tensor<double>& x) {
76 solve();
77 return 0.0;
78 }
79
80 void output_calc_info_schema(const std::string model, const double& energy) const;
81
82
83 std::string name() const {return "CC2";};
84
85 static void help() {
86 print_header2("help page for CC2");
87 print("The CC2 code computes correlated ground and excited state energies:\n");
88 print(" - MP2 ground state");
89 print(" - CC2 ground and excited states");
90 print(" - ADC(2) and CIS(D) excited states\n");
91 print("You need a SCF reference calculation from the nemo program. If there no such calculation can");
92 print("be found CC2 will perform its own. If excited states are requested also a CIS calculation is ");
93 print("necessary.\n");
94 print("Note that for correlated calculations the k parameter must be chosen small, typically k=5 or k=6 ");
95 print("because the curse of dimensions make higher k extremely expensive\n");
96 print("You can print all available calculation parameters by running\n");
97 print("cc2 --print_parameters\n");
98 print("You can perform a simple MP2 calculation by running\n");
99 print("cc2 --geometry=h2o.xyz\n");
100 print("provided you have an xyz file in your directory.\n");
101
102 }
103
104 static void print_parameters() {
106 print("\ndefault parameters for the cc2 program are\n");
107 param.print("cc2","end");
108 print("\n\nthe molecular geometry must be specified in a separate block:");
110 }
111
112 virtual bool selftest() {
113 return true;
114 };
115 void
116 plot(const real_function_3d& f, const std::string& msg = "unspecified function") const {
117 plot_plane(world, f, msg);
118 output("Plotted " + msg);
119 }
120
121 /// The World
123 /// Structure holds all the parameters used in the CC2 calculation
125 /// The SCF Calculation
126 std::shared_ptr<Nemo> nemo;
127 /// The excited state cis calculation
128 std::shared_ptr<TDHF> tdhf;
129 /// The CC Operator Class
130 CCPotentials CCOPS;
131 /// Formated Output (same as used in CC2Potentials structure)
133 /// map Pair struct to vector
135
136 /// solve the CC2 ground state equations, returns the correlation energy
137 void solve();
138
139
140 std::vector<CC_vecfunction>
141 solve_ccs() const;
142
143 double compute_mp3(const Pairs<CCPair>& mp2pairs) const {
144 MP3 mp3(CCOPS);
145 double mp3_contribution=mp3.mp3_energy_contribution_macrotask_driver(mp2pairs);
146 return mp3_contribution;
147 }
148
149 double
151
152 /// solve the excited state LR-CC2 equations for a given excitation
153
154 /// @param[in] gs_doubles: the ground state doubles
155 /// @param[in] gs_singles: the ground state singles
156 /// @param[in] cis: the CIS singles
157 /// @param[in] excitation: the excitation number
158 /// @return a tuple with the excited state doubles, the excited state singles and the excitation energy
159 std::tuple<Pairs<CCPair>, CC_vecfunction, double>
160 solve_lrcc2(Pairs<CCPair>& gs_doubles, const CC_vecfunction& gs_singles, const CC_vecfunction& cis,
161 const std::size_t excitation, Info& info) const;
162
163 double
165
166 /// convencience function to iterate the CC2 ground state singles,
167 /// makes the right call on the iterate_singles functions
168 static bool
170 // CCOPS.clear_potentials(singles);
172 Pairs<CCPair> empty;
174 empty, CT_CC2, info.parameters.iter_max_3D(), info);
175 }
176
177 bool
184
185 static bool
194
195 /// convenience function to iterate the CCS Response singles,
196 /// makes the right call on the iterate_singles functions
197 bool
199 Pairs<CCPair> empty;
200 // CCOPS.clear_potentials(x);
202 return iterate_singles(world, x, CC_vecfunction(PARTICLE), empty, empty, CT_LRCCS, info.parameters.iter_max_3D(), info);
203 }
204
205 /// Iterates the singles equations for CCS, CC2, LRCC2
206 /// The corresponding regulairzation tails of the doubles are updated in every iteration (therefore doubles are not marked as const)
207 /// @param[in] : singles, the singles that are iterated
208 /// @param[in] : singles2, additional set of singles for potential (in LRCC2 this are the Ground State singles)
209 /// @param[in] : gs_doubles, Ground State doubles (Needed for CC2 and LRCC2)
210 /// @param[in] : ex_doubles, Excited State doubles (Needed for LRCC2)
211 /// @param[in] : ctype: the calculation type: CCS, CC2, CC2_response_
212 /// @param[in] : maxiter: maxmial number of iterations
213 /// @param[out]: true if the overall change of the singles is below 10*donv_6D
214 static bool
216 const Pairs<CCPair>& ex_doubles, const CalcType ctype, const std::size_t maxiter, Info& info);
217
218 /// return the file name for singles
219 static std::string singles_name(const CalcType& ctype, const FuncType& type, int ex=-1) {
220 std::string fname=assign_name(ctype)+"_"+madness::name(type,ex);
221 return fname;
222 }
223
224 /// read singles from file or initialize new ones
225
226 /// type: PARTICLE (cc2) or RESPONSE (lrcc2)
227 /// default_to_zero: if true, initialize with zero functions, otherwise return empty vector
228 /// ex: if type is RESPONSE, the excitation number
229 static CC_vecfunction
230 initialize_singles(World&, const CalcType& ctype, const FuncType type, int ex=-1);
231
232 static CC_vecfunction
233 initialize_singles_to_zero(World& world, const CalcType& ctype, const FuncType type, const Info& info);
234
235 /// read pairs from file or initialize new ones
236 bool initialize_pairs(Pairs<CCPair>& pairs, const CCState ftype, const CalcType ctype, const CC_vecfunction& tau,
237 const CC_vecfunction& x, const size_t extitation, const Info& info) const;
238
239 /// Iterates a pair of the CC2 doubles equations
240 bool
242
243 bool
245
246 static bool
248 Pairs<CCPair>& lrcc2_d, const Info& info);
249
251 MADNESS_ASSERT(pair.ctype == CT_CISPD);
253 MADNESS_ASSERT(pair.bsh_eps == CCOPS.get_epsilon(pair.i, pair.j) + ccs.omega);
254 if (pair.constant_part.is_initialized()) return false; // the CIS(D) constant part does not change because there is no singles iteration (like MP2)
255 // make screening Operator
258 Gscreen.modified() = true;
259
260 if (parameters.QtAnsatz()) pair.constant_part = CCOPS.make_constant_part_cispd_Qt(pair, ccs, &Gscreen);
261 else pair.constant_part = CCOPS.make_constant_part_cispd(pair, ccs, &Gscreen);
262 save(pair.constant_part, pair.name() + "_const");
263 return true;
264
265 }
266
268 std::cout << assign_name(pair.ctype);
269 MADNESS_ASSERT(pair.ctype == CT_ADC2);
271 MADNESS_ASSERT(pair.bsh_eps == CCOPS.get_epsilon(pair.i, pair.j) + ccs.omega);
272 // make screening Operator
275 Gscreen.modified() = true;
276
277 if (parameters.QtAnsatz()) pair.constant_part = CCOPS.make_constant_part_cispd_Qt(pair, ccs, &Gscreen);
278 else pair.constant_part = CCOPS.make_constant_part_cispd(pair, ccs, &Gscreen);
279 save(pair.constant_part, pair.name() + "_const");
280 return true;
281
282 }
283
284 /// forward to the other function (converting CCPair to real_function)
285 static Pairs<real_function_6d> compute_local_coupling(const std::vector<CCPair> &vpairs, const Info& info) {
286 // create new pairs structure
287 Pairs<CCPair> pairs;
288 for (auto& tmp_pair : vpairs) pairs.insert(tmp_pair.i, tmp_pair.j, tmp_pair);
289 auto ccpair2function = [](const CCPair& a) {return a.function();};
291 };
292
293 /// compute the coupling of singles function if orbitals are localized
294
295 /// @return the coupling terms c_i = -\sum_(j\neq i) f_ij |\phi_j> (for whatever phi is)
296 static std::vector<real_function_3d> compute_local_coupling(const std::vector<real_function_3d>& singles,
297 const Info& info) {
298
299 MADNESS_CHECK_THROW(singles.size()>0,"compute_local_coupling: singles vector is empty");
300 World& world=singles.front().world();
301 auto active=Slice(info.parameters.freeze(),-1);
302 Tensor<double> Fact=info.fock(active,active);
303 for (int i=0; i<Fact.dim(0); ++i) Fact(i,i)=0.0;
305 return fock_coupling;
306 }
307
308 /// add the coupling terms for local MP2
309
310 /// \sum_{k\neq i} f_ki |u_kj> + \sum_{l\neq j} f_lj |u_il>
312
313
315
317
318 /// @return the new fock matrix
320};
321
322
323} /* namespace madness */
324
325#endif /* CC2_H_ */
Definition CC2.h:28
double solve_cispd(Pairs< CCPair > &doubles, const Pairs< CCPair > &mp2_pairs, const CC_vecfunction &cis_singles)
Definition CC2.cc:616
CCMessenger & output
Formated Output (same as used in CC2Potentials structure)
Definition CC2.h:132
static CC_vecfunction initialize_singles_to_zero(World &world, const CalcType &ctype, const FuncType type, const Info &info)
Definition CC2.cc:1328
Tensor< double > enforce_core_valence_separation(const Tensor< double > &fmat)
Definition CC2.cc:316
bool iterate_ccs_singles(CC_vecfunction &x, Info &info) const
Definition CC2.h:198
static CC_vecfunction initialize_singles(World &, const CalcType &ctype, const FuncType type, int ex=-1)
read singles from file or initialize new ones
Definition CC2.cc:1311
static std::string singles_name(const CalcType &ctype, const FuncType &type, int ex=-1)
return the file name for singles
Definition CC2.h:219
virtual bool selftest()
Definition CC2.h:112
double solve_mp2_coupled(Pairs< CCPair > &doubles, Info &info)
Definition CC2.cc:376
void output_calc_info_schema(const std::string model, const double &energy) const
Definition CC2.cc:296
World & world
The World.
Definition CC2.h:122
static void help()
Definition CC2.h:85
CCParameters parameters
Structure holds all the parameters used in the CC2 calculation.
Definition CC2.h:124
double compute_mp3(const Pairs< CCPair > &mp2pairs) const
Definition CC2.h:143
void solve()
solve the CC2 ground state equations, returns the correlation energy
Definition CC2.cc:19
static bool iterate_cc2_singles(World &world, CC_vecfunction &singles, Pairs< CCPair > &doubles, Info &info)
Definition CC2.h:169
std::tuple< Pairs< CCPair >, CC_vecfunction, double > solve_lrcc2(Pairs< CCPair > &gs_doubles, const CC_vecfunction &gs_singles, const CC_vecfunction &cis, const std::size_t excitation, Info &info) const
solve the excited state LR-CC2 equations for a given excitation
Definition CC2.cc:890
bool update_constant_part_adc2(const CC_vecfunction &ccs, CCPair &pair)
Definition CC2.h:267
bool update_constant_part_cispd(const CC_vecfunction &ccs, CCPair &pair)
Definition CC2.h:250
CC2(World &world_, const commandlineparser &parser, const std::shared_ptr< Nemo > nemo_)
Definition CC2.h:31
bool iterate_adc2_singles(Pairs< CCPair > &mp2, CC_vecfunction &singles, Pairs< CCPair > &x, Info &info)
Definition CC2.h:178
static void print_parameters()
Definition CC2.h:104
static Pairs< real_function_6d > compute_local_coupling(const std::vector< CCPair > &vpairs, const Info &info)
forward to the other function (converting CCPair to real_function)
Definition CC2.h:285
std::shared_ptr< Nemo > nemo
The SCF Calculation.
Definition CC2.h:126
double value()
Definition CC2.h:71
double solve_cc2(CC_vecfunction &tau, Pairs< CCPair > &u, Info &info) const
Definition CC2.cc:763
std::vector< CC_vecfunction > solve_ccs() const
Definition CC2.cc:361
virtual ~CC2()
Definition CC2.h:68
PairVectorMap triangular_map
map Pair struct to vector
Definition CC2.h:134
static bool iterate_singles(World &world, CC_vecfunction &singles, const CC_vecfunction singles2, const Pairs< CCPair > &gs_doubles, const Pairs< CCPair > &ex_doubles, const CalcType ctype, const std::size_t maxiter, Info &info)
Definition CC2.cc:1084
static bool iterate_lrcc2_pairs(World &world, const CC_vecfunction &cc2_s, const CC_vecfunction lrcc2_s, Pairs< CCPair > &lrcc2_d, const Info &info)
Definition CC2.cc:671
static std::vector< real_function_3d > compute_local_coupling(const std::vector< real_function_3d > &singles, const Info &info)
compute the coupling of singles function if orbitals are localized
Definition CC2.h:296
CCPotentials CCOPS
The CC Operator Class.
Definition CC2.h:130
bool check_core_valence_separation(const Tensor< double > &fmat) const
Definition CC2.cc:308
static bool iterate_lrcc2_singles(World &world, const CC_vecfunction &cc2_s, Pairs< CCPair > &cc2_d, CC_vecfunction &lrcc2_s, Pairs< CCPair > lrcc2_d, Info &info)
Definition CC2.h:186
double value(const Tensor< double > &x)
Should return the value of the objective function.
Definition CC2.h:75
bool iterate_pair(CCPair &pair, const CC_vecfunction &singles=CC_vecfunction(UNDEFINED)) const
Iterates a pair of the CC2 doubles equations.
Definition CC2.cc:977
bool iterate_adc2_pairs(Pairs< CCPair > &cispd, const CC_vecfunction &ccs)
Definition CC2.cc:654
std::shared_ptr< TDHF > tdhf
The excited state cis calculation.
Definition CC2.h:128
std::string name() const
Definition CC2.h:83
void plot(const real_function_3d &f, const std::string &msg="unspecified function") const
Definition CC2.h:116
bool initialize_pairs(Pairs< CCPair > &pairs, const CCState ftype, const CalcType ctype, const CC_vecfunction &tau, const CC_vecfunction &x, const size_t extitation, const Info &info) const
read pairs from file or initialize new ones
Definition CC2.cc:1340
Definition CCStructures.h:1180
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:100
static void set_truncate_mode(int value)
Sets the default truncation mode.
Definition funcdefaults.h:235
static void set_thresh(double value)
Sets the default threshold.
Definition funcdefaults.h:183
static void set_special_level(int value)
Existing functions are unaffected.
Definition funcdefaults.h:206
static void set_k(int value)
Sets the default wavelet order.
Definition funcdefaults.h:170
static void set_initial_level(int value)
Sets the default initial projection level.
Definition funcdefaults.h:200
Definition mp3.h:21
static void print_parameters()
Definition molecule.cc:110
class implementing properties of QC models
Definition QCPropertyInterface.h:11
Convolutions in separated form (including Gaussian)
Definition operator.h:139
A slice defines a sub-range or patch of a dimension.
Definition slice.h:103
Definition TDHF.h:24
A tensor is a multidimensional array.
Definition tensor.h:317
A parallel world class.
Definition world.h:132
double(* energy)()
Definition derivatives.cc:58
const int maxiter
Definition gygi_soltion.cc:68
static double u(double r, double c)
Definition he.cc:20
Implements (2nd generation) static load/data balancing for functions.
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:207
Main include file for MADNESS and defines Function interface.
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
CalcType
Calculation Types used by CC2.
Definition CCStructures.h:27
@ CT_CC2
Definition CCStructures.h:28
@ CT_LRCC2
Definition CCStructures.h:28
@ CT_LRCCS
Definition CCStructures.h:28
@ CT_CISPD
Definition CCStructures.h:28
@ CT_ADC2
Definition CCStructures.h:28
void print_header2(const std::string &s)
medium section heading
Definition print.cc:54
CCState
Type of Pairs used by CC_Pair2 class.
Definition CCStructures.h:31
@ EXCITED_STATE
Definition CCStructures.h:32
FuncType
Definition ccpairfunction.h:26
@ RESPONSE
Definition ccpairfunction.h:26
@ UNDEFINED
Definition ccpairfunction.h:26
@ PARTICLE
Definition ccpairfunction.h:26
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition vmra.h:689
std::vector< real_function_3d > vector_real_function_3d
Definition functypedefs.h:94
void plot_plane(World &world, const Function< double, NDIM > &function, const std::string name)
Definition funcplot.h:621
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition print.h:225
NDIM & f
Definition mra.h:2448
std::string type(const PairType &n)
Definition PNOParameters.h:18
std::string assign_name(const CCState &input)
Assigns strings to enums for formated output.
Definition CCStructures.cc:379
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
std::string name(const FuncType &type, const int ex=-1)
Definition ccpairfunction.h:28
void save(const Function< T, NDIM > &f, const std::string name)
Definition mra.h:2777
static const double a
Definition nonlinschro.cc:118
Implementation of Krylov-subspace nonlinear equation solver.
Implements most functionality of separated operators.
void clear_response()
clears only potentials of the response
Definition CCStructures.h:1026
void clear_all()
deltes all stored potentials
Definition CCStructures.h:1016
Definition CCStructures.h:90
void section(const std::string &msg) const
Definition CCStructures.cc:23
Definition CCStructures.h:213
std::size_t iter_max_3D() const
Definition CCStructures.h:329
double thresh_bsh_6D() const
Definition CCStructures.h:311
double thresh_3D() const
Definition CCStructures.h:301
double lo() const
Definition CCStructures.h:297
long freeze() const
Definition CCStructures.h:369
double thresh_6D() const
Definition CCStructures.h:305
bool QtAnsatz() const
Definition CCStructures.h:363
double dmin() const
Definition CCStructures.h:299
void information(World &world) const
print out the parameters
Definition CCStructures.cc:191
void sanity_check(World &world) const
check if parameters are set correct
Definition CCStructures.cc:202
A helper structure which holds a map of functions.
Definition CCStructures.h:523
POD holding some basic functions and some intermediates for the CC2 calculation.
Definition CCStructures.h:1108
Tensor< double > fock
Definition CCStructures.h:1114
CCIntermediatePotentials intermediate_potentials
Definition CCStructures.h:1115
CCParameters parameters
Definition CCStructures.h:1112
The interface to be provided by functions to be optimized.
Definition solvers.h:176
Definition CCStructures.h:395
Definition CCStructures.h:435
Pairs< R > convert(const Pairs< T > arg, const opT op) const
convert Pairs<T> to another type
Definition CCStructures.h:444
void insert(int i, int j, const T &pair)
Definition CCStructures.h:485
very simple command line parser
Definition commandlineparser.h:15
InputParameters param
Definition tdse.cc:203
Defines operations on vectors of Functions.