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 
18 namespace 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;
26  double log_high, log_low;
27  bool square_denominator=false;
28 
29  divide_add_interpolate(double hi, double lo, double eps_regularize) :
31  log_high(log10(thresh_high)), log_low(log10(thresh_low)) {}
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
53  ITERATOR(
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)
80 public:
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"));
114  if (dens_thresh_hi()<dens_thresh_lo()) {
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 
144 class OEP : public Nemo {
145 
146 private:
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 
157 public:
158 
159  OEP(World& world, const commandlineparser& parser)
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 
266  void load_restartdata(Tensor<double>& fock);
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);
326  real_function_3d lra=-1.0/(param.nalpha()+param.nbeta())*J.compute_potential(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) {
361  real_derivative_3d D(world,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
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
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
NemoCalculationParameters param
Definition: nemo.h:571
std::shared_ptr< SCF > get_calc() const
Definition: nemo.h:482
std::shared_ptr< SCF > calc
Definition: nemo.h:568
functor for a local U1 dot U1 potential
Definition: correlationfactor.h:532
Definition: oep.h:79
OEP_Parameters(World &world, const commandlineparser &parser)
Definition: oep.h:104
std::vector< double > kain_param() const
Definition: oep.h:139
unsigned int maxiter() const
Definition: oep.h:130
double dens_thresh_inv() const
Definition: oep.h:135
std::vector< std::string > model() const
Definition: oep.h:126
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
double levelshift() const
Definition: oep.h:131
double dens_thresh_lo() const
Definition: oep.h:134
Definition: oep.h:144
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
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
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
std::shared_ptr< Nemo > get_reference() const
Definition: oep.h:203
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::vector< double > compute_energy(const vecfuncT &nemo, const double E_X) const
Definition: oep.h:584
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
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
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
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
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:1436
void print_header2(const std::string &s)
medium section heading
Definition: print.cc:54
double abs(double x)
Definition: complexfun.h:48
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 truncate(World &world, response_space &v, double tol, bool fence)
Definition: basic_operators.cc:30
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
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
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
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< T, NDIM > sum(World &world, const std::vector< Function< T, NDIM > > &f, bool fence=true)
Returns new function — q = sum_i f[i].
Definition: vmra.h:1421
std::vector< Function< T, NDIM > > grad(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
shorthand gradient operator
Definition: vmra.h:1818
std::vector< Function< T, NDIM > > grad_bspline_one(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition: vmra.h:1878
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 > > 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:1667
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
void save(const Function< T, NDIM > &f, const std::string name)
Definition: mra.h:2745
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
double log_low
Definition: oep.h:26
std::vector< Tensor< double > > operator()(const Key< 3 > &key, const std::vector< Tensor< double > > &t) const
Definition: oep.h:39
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