MADNESS 0.10.1
oep.h
Go to the documentation of this file.
1/*
2 * oep.h
3 *
4 * Created on: Nov 6, 2019
5 * Author: fbischoff
6 */
7
8#ifndef SRC_APPS_CHEM_OEP_H_
9#define SRC_APPS_CHEM_OEP_H_
10
11
12
13#include<madness/chem/nemo.h>
17
18namespace madness {
19
20/// Class to compute terms of the potential
22
23 double thresh_high=1.e-5;
24 double thresh_low=1.e-7;
25 double eps_regularize=1.e-8;
28
32
33 std::size_t get_result_size() const {
34 return 1;
35 }
36
37 /// @param[in] t vector containing: oep, refdensity, longrange
38 /// @return (num1/denom1 - num2/denom2) * mask + (1-mask)*longrange
39 std::vector<Tensor<double> > operator()(const Key<3> & key,
40 const std::vector<Tensor<double> >& t) const {
41
42 const Tensor<double>& refdens=t[0];
43 const Tensor<double>& num1=t[1];
44 const Tensor<double>& denom1=t[2];
45 const Tensor<double>& num2=t[3];
46 const Tensor<double>& denom2=t[4];
47 const Tensor<double>& longrange=t[5];
48
49 // set 1 if above munge interval (dens > thresh_high) or 0 if below munge interval (dens < thresh_low)
50 // interpolate inbetween with logarithmic form of (r - lo)/(hi - lo) because this yields near linear behavior
51 Tensor<double> U = copy(refdens); // copy refdens in oder to have the same dimension in the tensor
52 U.fill(1.0); // start with a fuction that is 1.0 everywhere
54 U,
55 double r = refdens(IND);
56 double result=num1(IND)/(denom1(IND)+eps_regularize)
57 - num2(IND)/(denom2(IND)+eps_regularize);
58 if (square_denominator) result=num1(IND)/(denom1(IND)*denom1(IND)+eps_regularize)
59 - num2(IND)/(denom2(IND)*denom2(IND)+eps_regularize);
60 if (r > thresh_high) {
61 U(IND) = result;
62 } else if (r < thresh_low) {
63 U(IND) = longrange(IND);
64 } else {
65 // mask = 1 if refdens>hi, 0 if refdens < lo, interpolate in between
66 double mask=(log10(refdens(IND)) - log_low)/(log_high - log_low);
67 U(IND)=mask*result + (1.0-mask)*longrange(IND);
68 }
69 );
70 return std::vector<Tensor<double> > (1,U);
71 }
72};
73
74/// TODO:
75/// - change algorithm to handle localized orbitals
76/// - do not recompute the OEP potentials that don't depend on the orbital energies/Fock matrix elements
77/// - do not recompute the HF contributions to the OEP potentials
78/// - think about the long-range part of the Slater potential (or medium-range)
80public:
82 initialize<std::vector<std::string> >("model",{"dcep"},"comment on this: oaep ocep dcep mrks");
83 initialize<unsigned int>("maxiter",150,"maximum number of iterations in OEP algorithm");
84 initialize<bool>("restart",false,"restart from previous OEP calculation");
85 initialize<bool>("no_compute",false,"read from previous OEP calculation, no computation");
86 initialize<double>("levelshift",0.0,"shift occupied orbital energies in the BSH operator");
87// initialize<double>("conv_threshold",1.e-5,"comment on this");
88 initialize<double>("density_threshold_high",1.e-6,"comment on this");
89 initialize<double>("density_threshold_low",1.e-8,"comment on this");
90 initialize<double>("density_threshold_inv",1.e-9,"comment on this");
91 initialize<std::vector<double> >("kain_param",{1.0e-8, 3.0},"comment on this");
92
93// std::vector<bool> oep_model = {false, false, false, false};
94 initialize<unsigned int>("saving_amount",0,"choose level 0, 1, 2 or 3 for saving functions");
95 initialize<unsigned int>("save_iter_orbs",0,"if > 0 save all orbitals every ... iterations (needs a lot of storage!");
96 initialize<unsigned int>("save_iter_density",0,"if > 0 save KS density every ... iterations");
97 initialize<unsigned int>("save_iter_IKS",0,"if > 0 save IKS every ... iterations");
98 initialize<unsigned int>("save_iter_kin_tot_KS",0,"if > 0 save kin_tot_KS every ... iterations");
99 initialize<unsigned int>("save_iter_kin_P_KS",0,"if > 0 save kin_P_KS every ... iterations");
100 initialize<unsigned int>("save_iter_corrections",0,"if > 0 save OEP correction(s) every ... iterations");
101 initialize<unsigned int>("save_iter_effective_potential",0,"if > 0 save effective potential every ... iterations");
102 }
103
105 read_input_and_commandline_options(world, parser, "oep");
106 }
107
108 OEP_Parameters(const OEP_Parameters& other) = default;
109
110
112 set_derived_value("density_threshold_high",10.0*nparam.econv());
113 set_derived_value("density_threshold_low",0.01*get<double>("density_threshold_high"));
115 MADNESS_EXCEPTION("confused thresholds for the long-range transition",1);
116 }
117 set_derived_value("density_threshold_inv",0.1*get<double>("density_threshold_low"));
118
119 if (levelshift()>0.0) {
120 print("levelshift > 0.0 in oep parameters\n\n");
121 MADNESS_EXCEPTION("levelshift > 0.0 in oep parameters",1);
122 }
123 }
124
125 // convenience functions
126 std::vector<std::string> model() const {return get<std::vector<std::string> >("model");}
127
128 bool restart() const {return get<bool>("restart");}
129 bool no_compute() const {return get<bool>("no_compute");}
130 unsigned int maxiter() const {return get<unsigned int>("maxiter");}
131 double levelshift() const {return get<double>("levelshift");}
132// double conv_thresh() const {return get<double>("conv_threshold");}
133 double dens_thresh_hi() const {return get<double>("density_threshold_high");}
134 double dens_thresh_lo() const {return get<double>("density_threshold_low");}
135 double dens_thresh_inv() const{return get<double>("density_threshold_inv");}
136 unsigned int saving_amount() const {return get<unsigned int>("saving_amount");}
137 unsigned int save_iter_corrections () const {return get<unsigned int>("save_iter_corrections");}
138
139 std::vector<double> kain_param() const {return get<std::vector<double> >("kain_param");}
140
141};
142
143
144class OEP : public Nemo {
145
146private:
147
148 /// parameters for this OEP calculation
150
151 /// the wave function reference that determines the local potential
152 std::shared_ptr<Nemo> reference;
153
154 /// the final local potential
156
157public:
158
160 : Nemo(world, parser),
161 oep_param(world, parser) {
162
163 // add tight convergence criteria
164 std::vector<std::string> convergence_crit=param.get<std::vector<std::string> >("convergence_criteria");
165 if (std::find(convergence_crit.begin(),convergence_crit.end(),"each_energy")==convergence_crit.end()) {
166 convergence_crit.push_back("each_energy");
167 }
168 param.set_derived_value("convergence_criteria",convergence_crit);
169 calc->param.set_derived_value("convergence_criteria",convergence_crit);
170
171 // set reference
172 set_reference(std::make_shared<Nemo>(world,parser));
173 reference->param.set_derived_value("convergence_criteria",convergence_crit);
174 reference->get_calc()->param.set_derived_value("convergence_criteria",convergence_crit);
175
177 }
178
179 std::string name() const {return "oep";}
180
181 static void help() {
182 print_header2("help page for oep");
183 print("The oep code computes local exchange potentials based on a Hartree-Fock calculation from nemo");
184 print("oep --print_parameters\n");
185 print("You can perform a simple calculation by running\n");
186 print("oep --geometry=h2o.xyz\n");
187 print("provided you have an xyz file in your directory.");
188
189 }
190
191 static void print_parameters() {
193 print("default parameters for the oep program are");
194 param.print("oep", "end");
195 print("\n\nthe molecular geometry must be specified in a separate block:");
197 }
198
199 void set_reference(const std::shared_ptr<Nemo> reference1) {
200 reference=reference1;
201 }
202
203 std::shared_ptr<Nemo> get_reference() const {
204 return reference;
205 }
206
207 void print_parameters(std::vector<std::string> what) const {
208 for (auto w : what) {
209 if (w=="oep") oep_param.print("oep");
210 else if (w=="reference") reference->param.print("dft");
211 else if (w=="oep_calc") param.print("oep_calc");
212 else {MADNESS_EXCEPTION(std::string("unknown parameter set to print "+w).c_str(),1);}
213 }
214 }
215
217 return Vfinal;
218 }
219
220 virtual double value() {
221 return value(calc->molecule.get_all_coords());
222 }
223
224 /// update the json file with calculation input and output
225 void output_calc_info_schema(const double& energy) const;
226
227
228 virtual double value(const Tensor<double>& x) {
229 reference->value();
231 calc->copy_data(world,*(reference->get_calc()));
232 double energy=0.0;
233 Tensor<double> fock;
234 bool load_mos=(oep_param.restart() or oep_param.no_compute());
235 if (load_mos) load_restartdata(fock);
236
237 if (not oep_param.no_compute()) energy=solve(reference->get_calc()->get_amo());
239 return energy;
240 };
241
242 /// Iterative energy calculation for approximate OEP with EXACT EXCHANGE functional
243 /// for other functionals, slater potential must be modified
244 /// HF orbitals and eigenvalues are used as the guess here
245 /// note that KS_nemo is a reference and changes oep->get_calc()->amo orbitals
246 /// same for orbital energies (eigenvalues) KS_eigvals which is oep->get_calc()->aeps
247 /// converged if norm, total energy difference and orbital energy differences (if not OAEP) are converged
248 double solve(const vecfuncT& HF_nemo);
249
250 void analyze();
251
252 double iterate(const std::string model, const vecfuncT& HF_nemo, const tensorT& HF_eigvals,
253 vecfuncT& KS_nemo, tensorT& KS_Fock, real_function_3d& Voep,
254 const real_function_3d Vs) const;
255
256 virtual std::shared_ptr<Fock<double,3>> make_fock_operator() const;
257
258// MolecularOrbitals<double,3> to_MO() const {
259// std::vector<std::string> str_irreps;
260// vecfuncT aaa=symmetry_projector(calc->amo,R_square,str_irreps);
261// return MolecularOrbitals<double,3>(aaa,this->get_calc()->aeps,str_irreps,this->get_calc()->aocc,this->get_calc()->aset);
262// }
263
264 void save_restartdata(const Tensor<double>& fock) const;
265
267
268 std::tuple<Tensor<double>, vecfuncT> recompute_HF(const vecfuncT& HF_nemo) const;
269
270 /// The following function tests all essential parts of the OEP program qualitatively and some also quantitatively
271 bool selftest();
272
273 bool need_ocep_correction(const std::string& model) const {
274 return (model=="ocep") or (model=="dcep") or (model=="mrks");
275 }
276
277 bool need_dcep_correction(const std::string& model) const {
278 return (model=="dcep");
279 }
280
281 bool need_mrks_correction(const std::string& model) const {
282 return (model=="mrks");
283 }
284
285 /// print orbital energies in reverse order with optional shift
286 void print_orbens(const tensorT orbens) const {
287 for (long i = orbens.size() - 1; i >= 0; i--) {
288 printf(" e%2.2lu = %12.8f Eh\n", i, orbens(i));
289 }
290 }
291
292 double compute_and_print_final_energies(const std::string model, const real_function_3d& Voep,
293 const vecfuncT& KS_nemo, const tensorT& KS_Fock,
294 const vecfuncT& HF_nemo, const tensorT& HF_Fock) const;
295
296 /// compute density from orbitals with ragularization (Bischoff, 2014_1, equation (19))
298 real_function_3d density = 2.0*R_square*dot(world, nemo, nemo); // 2 because closed shell
299 return density;
300 }
301
302 /// compute Delta rho as an indicator for the result's quality
303 double compute_delta_rho(const real_function_3d rho_HF, const real_function_3d rho_KS) const {
304 // from Ospadov_2017, equation (26)
305 real_function_3d rho_diff = abs(rho_KS - rho_HF);
306 double Drho = rho_diff.trace();
307 return Drho;
308 }
309
310 /// compute Slater potential (Kohut, 2014, equation (15))
312
313 Exchange<double,3> K(world,reference->get_calc()->param.lo());
314 K.set_bra_and_ket(R_square * nemo, nemo);
315 const vecfuncT Knemo = K(nemo);
316 // 2.0*R_square in numerator and density (rho) cancel out upon division
317 real_function_3d numerator = -1.0*dot(world, nemo, Knemo);
318 real_function_3d rho = dot(world, nemo, nemo);
319
320 // long-range asymptotic behavior for Slater potential is \int 1/|r-r'| * |phi_HOMO|^2 dr'
321 // in order to compute this lra, use Coulomb potential with only HOMO density (= |phi_HOMO|^2)
322// Coulomb J(world);
323// J.reset_poisson_operator_ptr(param.lo(),param.econv());
324// real_function_3d lra = -1.0*J.compute_potential(R_square*square(nemo[homo_ind]));
325 Coulomb<double,3> J(world,this);
327// print("compute long-range part of the Slater potential from the full molecular density");
328 if (oep_param.saving_amount() >= 3) save(lra, "lra_slater");
329
331 real_function_3d one=real_factory_3d(world).functor([](const coord_3d& r){return 1.0;});
332 std::vector<real_function_3d> args={R_square*rho,numerator,rho,zero,one,lra};
334
338 return VSlater;
339 }
340
341
342 /// compute the total kinetic energy density of equation (6) from Kohut
343
344 /// return without the NCF and factor 2 for closed shell !
346
347 // compute tau = 1/2 * sum |\nabla phi_i|^2
348 // = 1/2 * sum {(\nabla R)^2 * nemo_i^2 + 2 * R * nemo_i * (\nabla R) * (\nabla nemo_i)) + R^2 * (\nabla nemo_i)^2}
349 // = 1/2 * R^2 * sum {U1dot * nemo_i^2 + 2 * nemo_i * U1 * (\nabla nemo_i)) + (\nabla nemo_i)^2}
350
351 // get \nabla R and (\nabla R)^2 via and U1 = -1/R \nabla R and U1dot = (1/R \nabla R)^2
352 const vecfuncT U1 = this->ncf->U1vec();
354 const real_function_3d U1dot = real_factory_3d(world).functor(u1_dot_u1).truncate_on_project();
355
356 real_function_3d u1u1term=U1dot*dot(world,nemo,nemo);
357 real_function_3d u1dnterm=real_factory_3d(world).compressed();
358 real_function_3d dndnterm=real_factory_3d(world).compressed();
359
360 for (int idim=0; idim<3; ++idim) {
362 if(param.dft_deriv() == "bspline") D.set_bspline1();
363 vecfuncT nemo_copy=copy(world,nemo);
364 refine(world,nemo_copy);
365 std::vector<real_function_3d> dnemo=apply(world,D,nemo_copy);
366 u1dnterm+=U1[idim]*dot(world,nemo_copy,dnemo);
367 refine(world,dnemo);
368 dndnterm+=dot(world,dnemo,dnemo);
369 }
370
371 real_function_3d tau=0.5*(u1u1term - 2.0*u1dnterm + dndnterm);
372
373 return tau;
374
375 }
376
377 /// compute the Pauli kinetic energy density divided by the density tau_P/rho with equation (16) from Ospadov, 2017
379
380 // compute the numerator tau_P
381
382 // get \nabla nemo
383 std::vector<vecfuncT> grad_nemo(nemo.size());
384 for (size_t i = 0; i < nemo.size(); i++) {
385 vecfuncT nemo_copy=copy(world,nemo);
386 refine(world,nemo_copy);
387
388 if(param.dft_deriv() == "bspline") grad_nemo[i] = grad_bspline_one(nemo_copy[i]); // gradient using b-spline
389 else grad_nemo[i] = grad(nemo_copy[i]); // default gradient using abgv
390 }
391
392 vecfuncT grad_nemo_term;
393 for (size_t i = 0; i < nemo.size(); i++) {
394 for (size_t j = i + 1; j < nemo.size(); j++) {
395 vecfuncT tmp = nemo[i]*grad_nemo[j] - nemo[j]*grad_nemo[i];
396 grad_nemo_term.push_back(dot(world, tmp, tmp));
397 }
398 }
399 real_function_3d numerator = sum(world, grad_nemo_term); // numerator = tau_P * 2 * rho / R^4
400
401 return numerator;
402
403 }
404
405 /// compute correction of the given model
406
407 /// shift of the diagonal elements of the fock matrix results in a global shift
408 /// of the potential
409 /// \f[
410 /// \frac{\bar \epsilon}{\rho} = \frac{1}{\rho}\sum_{ij}\phi_i(F_ij+\delta_ij s)\phi_j
411 /// = \frac{1}{\rho} ( \sum_{ij}\phi_i F_ij \phi_j + s\sum_i\phi_i\phi_i )
412 /// = s + \frac{1}{\rho} \sum_{ij}\phi_i F_ij \phi_j
413 /// \f]
415 const vecfuncT& nemoHF, const vecfuncT& nemoKS,
416 const tensorT& fockHF, const tensorT& fockKS) const {
417
418 if (fockKS.normf()<1.e-10) {
420 return zero;
421 }
422
423 // compute HOMO energies and the long-range asymptotics
424 auto [eval, evec] = syev(fockKS);
425 double homoKS = -eval.max();
426
427 auto [eval1, evec1] = syev(fockHF);
428 double homoHF = -eval1.max();
429
430 double longrange=homoHF-homoKS;
431
432 tensorT fock1=copy(fockKS);
433 for (int i=0; i<fock1.dim(0); ++i) fock1(i,i)-=longrange;
434
435 // 2.0*R_square in numerator and density (rho) cancel out upon division
436// real_function_3d numeratorHF=-1.0*compute_energy_weighted_density_local(nemoHF,fockHF);
437// real_function_3d numeratorHF=-1.0*compute_energy_weighted_density(nemoHF,eigvalsHF);
438 real_function_3d numeratorKS=-1.0*compute_energy_weighted_density_local(nemoKS,fock1);
439// real_function_3d numeratorKS=-1.0*compute_energy_weighted_density(nemoKS,eval);
440
441 // 2.0*R_square in numerator and density (rho) cancel out upon division
442 real_function_3d densityKS = dot(world, nemoKS, nemoKS);
443 real_function_3d densityHF = dot(world, nemoHF, nemoHF);
444
446 std::vector<real_function_3d> args={densityKS,ocep_numerator_HF,densityHF,
447 numeratorKS,densityKS,lra_func};
449
452 real_function_3d correction=multi_to_multi_op_values(op,args)[0];
453
454
455 return correction;
456
457 }
458
459 /// return without the NCF and factor 2 for closed shell !
460
461 /// follow Eq. (4) of Ryabinkin2014
463 const tensorT& fock) const {
464 return dot(world,nemo,transform(world,nemo,fock));
465 }
466
467 /// compute correction of the given model
469 const vecfuncT& nemoHF, const vecfuncT& nemoKS) const {
470
471 // 2.0*R_square in numerator and density (rho) cancel out upon division
472// real_function_3d numeratorHF=compute_total_kinetic_density(nemoHF);
474
475 // 2.0*R_square in numerator and density (rho) cancel out upon division
476 real_function_3d densityKS = dot(world, nemoKS, nemoKS);
477 real_function_3d densityHF = dot(world, nemoHF, nemoHF);
478
479 real_function_3d lra_func=real_factory_3d(world).functor([](const coord_3d& r)
480 {return 0.0;});
481
482 std::vector<real_function_3d> args={densityKS,dcep_numerator_HF,densityHF,
483 numeratorKS,densityKS,lra_func};
485
488 real_function_3d correction=multi_to_multi_op_values(op,args)[0];
489
490 return correction;
491 }
492
493 /// compute correction of the given model
495 const vecfuncT& nemoHF, const vecfuncT& nemoKS) const {
496
497 // 2.0*R_square in numerator and density (rho) cancel out upon division
498// real_function_3d numeratorHF=compute_Pauli_kinetic_density(nemoHF);
500
501 // 2.0*R_square in numerator and density (rho) cancel out upon division
502 real_function_3d densityKS = dot(world, nemoKS, nemoKS);
503 real_function_3d densityHF = dot(world, nemoHF, nemoHF);
504
505 // longrange correction is zero
506 real_function_3d lra_func=real_factory_3d(world).functor([](const coord_3d& r) {return 0.0;});
507
508 real_function_3d denssq=square(densityKS);
509 std::vector<real_function_3d> args={densityKS,mrks_numerator_HF,densityHF,
510 numeratorKS,densityKS,lra_func};
512
515 op.square_denominator=true;
516 real_function_3d correction=0.5*multi_to_multi_op_values(op,args)[0];
517
518 return correction;
519 }
520
521 /// compute all potentials from given nemos except kinetic energy
522 void compute_nemo_potentials(const vecfuncT& nemo, vecfuncT& Jnemo, vecfuncT& Unemo) const {
523
524 // compute Coulomb part
525 compute_coulomb_potential(nemo, Jnemo);
526
527 // compute nuclear potential part
528 Nuclear<double,3> Unuc(world, this->ncf);
529 Unemo = Unuc(nemo);
530
531 }
532
533 /// compute Coulomb potential
534 void compute_coulomb_potential(const vecfuncT& nemo, vecfuncT& Jnemo) const {
535
536 Coulomb<double,3> J(world, this);
539 Jnemo = J(nemo);
540 truncate(world, Jnemo);
541
542 }
543
544 /// compute exchange potential (needed for Econv)
545 void compute_exchange_potential(const vecfuncT& nemo, vecfuncT& Knemo) const {
546
548 K.set_bra_and_ket(R_square * nemo, nemo);
549 Knemo = K(nemo);
550 truncate(world, Knemo);
551
552 }
553
554 /// compute Evir using Levy-Perdew virial relation (Kohut_2014, (43) or Ospadov_2017, (25))
555 double compute_exchange_energy_vir(const vecfuncT& nemo, const real_function_3d Vx) const {
556
557 // make vector of functions r = (x, y, z)
558 auto monomial_x = [] (const coord_3d& r) {return r[0];};
559 auto monomial_y = [] (const coord_3d& r) {return r[1];};
560 auto monomial_z = [] (const coord_3d& r) {return r[2];};
561 vecfuncT r(3);
562 r[0]=real_factory_3d(world).functor(monomial_x);
563 r[1]=real_factory_3d(world).functor(monomial_y);
564 r[2]=real_factory_3d(world).functor(monomial_z);
565
566 real_function_3d rhonemo = 2.0*dot(world, nemo, nemo); // 2 because closed shell
567
568 real_function_3d rhoterm = 3*rhonemo + dot(world, r, -2.0*ncf->U1vec()*rhonemo + grad(rhonemo,true));
569 double Ex = inner(Vx, R_square*rhoterm);
570 return Ex;
571
572 }
573
574 /// compute exchange energy using the expectation value of the exchange operator
575 double compute_exchange_energy_conv(const vecfuncT phi, const vecfuncT Kphi) const {
576
577 double Ex = -1.0*inner(world, phi, Kphi).sum();
578 return Ex;
579
580 }
581
582 /// compute all energy contributions except exchange and sum up for total energy
583 /// the exchange energy must be computed priorly with the compute_exchange_energy_... methods
584 std::vector<double> compute_energy(const vecfuncT& nemo, const double E_X) const {
585
586 // compute kinetic energy
587 double E_kin = compute_kinetic_energy(nemo);
588
590
591 // compute external potential (nuclear attraction)
592 real_function_3d Vext = calc->potentialmanager->vnuclear();
593 Coulomb<double,3> J(world,this);
595
596 // compute remaining energies: nuclear attraction, Coulomb, nuclear repulsion
597 // computed as expectation values (see Szabo, Ostlund (3.81))
598 const double E_ext = inner(Vext,density); // 2 because closed shell
599 const double E_J = 0.5*inner(density,Jpotential);
600 const double E_nuc = calc->molecule.nuclear_repulsion_energy();
601 double energy = E_kin + E_ext + E_J + E_X + E_nuc;
602
603 if (world.rank() == 0) {
604 printf("\n kinetic energy %15.8f Eh\n", E_kin);
605 printf(" electron-nuclear attraction energy %15.8f Eh\n", E_ext);
606 printf(" Coulomb energy %15.8f Eh\n", E_J);
607 printf(" exchange energy %15.8f Eh\n", E_X);
608 printf(" nuclear-nuclear repulsion energy %15.8f Eh\n", E_nuc);
609 printf(" total energy %15.8f Eh\n\n", energy);
610 }
611 return {energy,E_kin,E_ext,E_J,E_X};
612
613 }
614
615 /// compute diagonal elements of Fock matrix
617 const vecfuncT& phi, const vecfuncT& Kphi, const real_function_3d& Vx) const {
618 return KS_eigvals - inner(world, phi, Kphi) - inner(world, phi, Vx*phi);
619 }
620
621 /// cumpute E^(0) = \sum_i \epsilon_i^KS
622 double compute_E_zeroth(const tensorT eigvals) const {
623 double E_0 = 0.0;
624 for (long i = 0; i < eigvals.size(); i++) {
625 E_0 += eigvals(i);
626 }
627 E_0 *= 2.0; // closed shell: every orbital energy must be counted twice
628 return E_0;
629 }
630
631 /// compute E^(1) = 1/2*\sum_ij <ij||ij> - \sum_i <i|J + Vx|i> = \sum_i <i|- 0.5*J - 0.5*K - Vx|i>
632 double compute_E_first(const vecfuncT phi, const vecfuncT Jphi, const vecfuncT Kphi, const real_function_3d Vx) const {
633
634 //compute expectation values:
635 const double E_J = inner(world, phi, Jphi).sum();
636 const double E_K = inner(world, phi, Kphi).sum();
637 const double E_Vx = inner(world, phi, Vx*phi).sum();
638 printf(" E_J = %15.8f Eh\n", E_J);
639 printf(" E_K = %15.8f Eh\n", E_K);
640 printf(" E_Vx = %15.8f Eh\n", E_Vx);
641 const double E_nuc = calc->molecule.nuclear_repulsion_energy();
642
643 double E_1 = -1.0*(E_J + E_K + 2.0*E_Vx) + E_nuc;
644 return E_1;
645
646 }
647
648};
649
650
651
652} /* namespace madness */
653
654#endif /* SRC_APPS_CHEM_OEP_H_ */
double w(double t, double eps)
Definition DKops.h:22
Operators for the molecular HF and DFT code.
long dim(int i) const
Returns the size of dimension i.
Definition basetensor.h:147
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
Definition SCFOperators.h:334
Function< T, NDIM > compute_potential(const Function< T, NDIM > &density) const
given a density compute the Coulomb potential
Definition SCFOperators.h:432
const real_function_3d & potential() const
getter for the Coulomb potential
Definition SCFOperators.h:421
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition derivative.h:266
Definition SCFOperators.h:104
Exchange & set_bra_and_ket(const vecfuncT &bra, const vecfuncT &ket)
Definition SCFOperators.cc:699
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition mra.h:1099
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:66
static void print_parameters()
Definition molecule.cc:110
real_function_3d R_square
the square of the nuclear correlation factor
Definition nemo.h:319
double compute_kinetic_energy(const std::vector< Function< T, NDIM > > &nemo) const
compute kinetic energy as square of the "analytical" expectation value
Definition nemo.h:183
World & world
Definition nemo.h:310
The Nemo class.
Definition nemo.h:326
void set_protocol(const double thresh)
adapt the thresholds consistently to a common value
Definition nemo.h:617
std::shared_ptr< SCF > get_calc() const
Definition nemo.h:482
NemoCalculationParameters param
Definition nemo.h:571
std::shared_ptr< SCF > calc
Definition nemo.h:568
functor for a local U1 dot U1 potential
Definition correlationfactor.h:532
Definition SCFOperators.h:455
Definition oep.h:79
OEP_Parameters(World &world, const commandlineparser &parser)
Definition oep.h:104
std::vector< std::string > model() const
Definition oep.h:126
unsigned int maxiter() const
Definition oep.h:130
double dens_thresh_inv() const
Definition oep.h:135
bool restart() const
Definition oep.h:128
OEP_Parameters(const OEP_Parameters &other)=default
bool no_compute() const
Definition oep.h:129
unsigned int save_iter_corrections() const
Definition oep.h:137
void set_derived_values(const Nemo::NemoCalculationParameters &nparam)
Definition oep.h:111
unsigned int saving_amount() const
Definition oep.h:136
double dens_thresh_hi() const
Definition oep.h:133
OEP_Parameters()
Definition oep.h:81
std::vector< double > kain_param() const
Definition oep.h:139
double levelshift() const
Definition oep.h:131
double dens_thresh_lo() const
Definition oep.h:134
Definition oep.h:144
bool selftest()
The following function tests all essential parts of the OEP program qualitatively and some also quant...
Definition madness/chem/oep.cc:402
bool need_mrks_correction(const std::string &model) const
Definition oep.h:281
void set_reference(const std::shared_ptr< Nemo > reference1)
Definition oep.h:199
OEP_Parameters oep_param
parameters for this OEP calculation
Definition oep.h:149
static void help()
Definition oep.h:181
double compute_exchange_energy_vir(const vecfuncT &nemo, const real_function_3d Vx) const
compute Evir using Levy-Perdew virial relation (Kohut_2014, (43) or Ospadov_2017, (25))
Definition oep.h:555
real_function_3d compute_total_kinetic_density(const vecfuncT &nemo) const
compute the total kinetic energy density of equation (6) from Kohut
Definition oep.h:345
void save_restartdata(const Tensor< double > &fock) const
Definition madness/chem/oep.cc:84
virtual std::shared_ptr< Fock< double, 3 > > make_fock_operator() const
the OEP Fock operator is the HF Fock operator without exchange but with the OEP
Definition madness/chem/oep.cc:132
Tensor< double > compute_fock_diagonal_elements(const Tensor< double > &KS_eigvals, const vecfuncT &phi, const vecfuncT &Kphi, const real_function_3d &Vx) const
compute diagonal elements of Fock matrix
Definition oep.h:616
double compute_and_print_final_energies(const std::string model, const real_function_3d &Voep, const vecfuncT &KS_nemo, const tensorT &KS_Fock, const vecfuncT &HF_nemo, const tensorT &HF_Fock) const
Definition madness/chem/oep.cc:140
double compute_exchange_energy_conv(const vecfuncT phi, const vecfuncT Kphi) const
compute exchange energy using the expectation value of the exchange operator
Definition oep.h:575
double compute_delta_rho(const real_function_3d rho_HF, const real_function_3d rho_KS) const
compute Delta rho as an indicator for the result's quality
Definition oep.h:303
void print_parameters(std::vector< std::string > what) const
Definition oep.h:207
void compute_nemo_potentials(const vecfuncT &nemo, vecfuncT &Jnemo, vecfuncT &Unemo) const
compute all potentials from given nemos except kinetic energy
Definition oep.h:522
real_function_3d compute_slater_potential(const vecfuncT &nemo) const
compute Slater potential (Kohut, 2014, equation (15))
Definition oep.h:311
virtual double value()
Definition oep.h:220
bool need_dcep_correction(const std::string &model) const
Definition oep.h:277
real_function_3d compute_energy_weighted_density_local(const vecfuncT &nemo, const tensorT &fock) const
return without the NCF and factor 2 for closed shell !
Definition oep.h:462
real_function_3d compute_mrks_correction(const real_function_3d &mrks_numerator_HF, const vecfuncT &nemoHF, const vecfuncT &nemoKS) const
compute correction of the given model
Definition oep.h:494
std::tuple< Tensor< double >, vecfuncT > recompute_HF(const vecfuncT &HF_nemo) const
Definition madness/chem/oep.cc:120
std::string name() const
Definition oep.h:179
virtual double value(const Tensor< double > &x)
Should return the value of the objective function.
Definition oep.h:228
std::shared_ptr< Nemo > get_reference() const
Definition oep.h:203
double compute_E_first(const vecfuncT phi, const vecfuncT Jphi, const vecfuncT Kphi, const real_function_3d Vx) const
compute E^(1) = 1/2*\sum_ij <ij||ij> - \sum_i <i|J + Vx|i> = \sum_i <i|- 0.5*J - 0....
Definition oep.h:632
std::shared_ptr< Nemo > reference
the wave function reference that determines the local potential
Definition oep.h:152
real_function_3d compute_Pauli_kinetic_density(const vecfuncT &nemo) const
compute the Pauli kinetic energy density divided by the density tau_P/rho with equation (16) from Osp...
Definition oep.h:378
double iterate(const std::string model, const vecfuncT &HF_nemo, const tensorT &HF_eigvals, vecfuncT &KS_nemo, tensorT &KS_Fock, real_function_3d &Voep, const real_function_3d Vs) const
Definition madness/chem/oep.cc:214
void compute_coulomb_potential(const vecfuncT &nemo, vecfuncT &Jnemo) const
compute Coulomb potential
Definition oep.h:534
void print_orbens(const tensorT orbens) const
print orbital energies in reverse order with optional shift
Definition oep.h:286
void load_restartdata(Tensor< double > &fock)
Definition madness/chem/oep.cc:102
double compute_E_zeroth(const tensorT eigvals) const
cumpute E^(0) = \sum_i \epsilon_i^KS
Definition oep.h:622
real_function_3d compute_density(const vecfuncT &nemo) const
compute density from orbitals with ragularization (Bischoff, 2014_1, equation (19))
Definition oep.h:297
static void print_parameters()
Definition oep.h:191
real_function_3d get_final_potential() const
Definition oep.h:216
real_function_3d Vfinal
the final local potential
Definition oep.h:155
OEP(World &world, const commandlineparser &parser)
Definition oep.h:159
void analyze()
Definition madness/chem/oep.cc:67
void compute_exchange_potential(const vecfuncT &nemo, vecfuncT &Knemo) const
compute exchange potential (needed for Econv)
Definition oep.h:545
void output_calc_info_schema(const double &energy) const
update the json file with calculation input and output
Definition madness/chem/oep.cc:58
std::vector< double > compute_energy(const vecfuncT &nemo, const double E_X) const
Definition oep.h:584
bool need_ocep_correction(const std::string &model) const
Definition oep.h:273
real_function_3d compute_dcep_correction(const real_function_3d &dcep_numerator_HF, const vecfuncT &nemoHF, const vecfuncT &nemoKS) const
compute correction of the given model
Definition oep.h:468
real_function_3d compute_ocep_correction(const real_function_3d &ocep_numerator_HF, const vecfuncT &nemoHF, const vecfuncT &nemoKS, const tensorT &fockHF, const tensorT &fockKS) const
compute correction of the given model
Definition oep.h:414
double solve(const vecfuncT &HF_nemo)
Definition madness/chem/oep.cc:22
class for holding the parameters for calculation
Definition QCCalculationParametersBase.h:290
virtual void read_input_and_commandline_options(World &world, const commandlineparser &parser, const std::string tag)
Definition QCCalculationParametersBase.h:325
T get(const std::string key) const
Definition QCCalculationParametersBase.h:299
void print(const std::string header="", const std::string footer="") const
print all parameters
Definition QCCalculationParametersBase.cc:22
void set_derived_value(const std::string &key, const T &value)
Definition QCCalculationParametersBase.h:403
virtual real_function_3d density() const
Definition QCPropertyInterface.h:25
A tensor is a multidimension array.
Definition tensor.h:317
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition tensor.h:1726
Tensor< T > & fill(T x)
Inplace fill with a scalar (legacy name)
Definition tensor.h:562
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:318
double(* energy)()
Definition derivatives.cc:58
static double lo
Definition dirac-hatom.cc:23
real_function_3d mask
Definition dirac-hatom.cc:27
vecfuncT K(vecfuncT &ket, vecfuncT &bra, vecfuncT &vf)
Tensor< double > op(const Tensor< double > &x)
Definition kain.cc:508
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition madness_exception.h:119
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
void print_header2(const std::string &s)
medium section heading
Definition print.cc:54
double abs(double x)
Definition complexfun.h:48
Function< T, NDIM > square(const Function< T, NDIM > &f, bool fence=true)
Create a new function that is the square of f - global comm only if not reconstructed.
Definition mra.h:2681
std::vector< Function< T, NDIM > > grad_bspline_one(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition vmra.h:1935
void truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
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< Function< T, NDIM > > multi_to_multi_op_values(const opT &op, const std::vector< Function< T, NDIM > > &vin, const bool fence=true)
apply op on the input vector yielding an output vector of functions
Definition vmra.h:1724
FunctionFactory< double, 3 > real_factory_3d
Definition functypedefs.h:93
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
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
void refine(World &world, const std::vector< Function< T, NDIM > > &vf, bool fence=true)
refine the functions according to the autorefine criteria
Definition vmra.h:208
Function< TENSOR_RESULT_TYPE(T, R), NDIM > dot(World &world, const std::vector< Function< T, NDIM > > &a, const std::vector< Function< R, NDIM > > &b, bool fence=true)
Multiplies and sums two vectors of functions r = \sum_i a[i] * b[i].
Definition vmra.h:1493
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
vector< functionT > vecfuncT
Definition corepotential.cc:58
void refine_to_common_level(World &world, std::vector< Function< T, NDIM > > &vf, bool fence=true)
refine all functions to a common (finest) level
Definition vmra.h:218
std::vector< Function< T, NDIM > > grad(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
shorthand gradient operator
Definition vmra.h:1875
void save(const Function< T, NDIM > &f, const std::string name)
Definition mra.h:2745
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition mra.h:2002
void syev(const Tensor< T > &A, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Real-symmetric or complex-Hermitian eigenproblem.
Definition lapack.cc:969
Definition test_ar.cc:204
double econv() const
Definition CalculationParameters.h:141
int nalpha() const
Definition CalculationParameters.h:162
std::string dft_deriv() const
Definition CalculationParameters.h:194
int nbeta() const
Definition CalculationParameters.h:163
double lo() const
Definition CalculationParameters.h:177
class holding parameters for a nemo calculation beyond the standard dft parameters from moldft
Definition nemo.h:333
very simple command line parser
Definition commandlineparser.h:15
Class to compute terms of the potential.
Definition oep.h:21
double log_high
Definition oep.h:26
std::vector< Tensor< double > > operator()(const Key< 3 > &key, const std::vector< Tensor< double > > &t) const
Definition oep.h:39
double log_low
Definition oep.h:26
bool square_denominator
Definition oep.h:27
std::size_t get_result_size() const
Definition oep.h:33
double eps_regularize
Definition oep.h:25
divide_add_interpolate(double hi, double lo, double eps_regularize)
Definition oep.h:29
double thresh_high
Definition oep.h:23
double thresh_low
Definition oep.h:24
Definition dirac-hatom.cc:108
#define ITERATOR(t, exp)
Definition tensor_macros.h:249
#define IND
Definition tensor_macros.h:204
AtomicInt sum
Definition test_atomicint.cc:46