8 #ifndef SRC_APPS_CHEM_OEP_H_
9 #define SRC_APPS_CHEM_OEP_H_
55 double r = refdens(
IND);
70 return std::vector<Tensor<double> > (1,U);
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");
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");
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");
113 set_derived_value(
"density_threshold_low",0.01*get<double>(
"density_threshold_high"));
117 set_derived_value(
"density_threshold_inv",0.1*get<double>(
"density_threshold_low"));
120 print(
"levelshift > 0.0 in oep parameters\n\n");
126 std::vector<std::string>
model()
const {
return get<std::vector<std::string> >(
"model");}
128 bool restart()
const {
return get<bool>(
"restart");}
130 unsigned int maxiter()
const {
return get<unsigned int>(
"maxiter");}
131 double levelshift()
const {
return get<double>(
"levelshift");}
136 unsigned int saving_amount()
const {
return get<unsigned int>(
"saving_amount");}
139 std::vector<double>
kain_param()
const {
return get<std::vector<double> >(
"kain_param");}
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");
169 calc->param.set_derived_value(
"convergence_criteria",convergence_crit);
173 reference->param.set_derived_value(
"convergence_criteria",convergence_crit);
174 reference->get_calc()->param.set_derived_value(
"convergence_criteria",convergence_crit);
179 std::string
name()
const {
return "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.");
193 print(
"default parameters for the oep program are");
195 print(
"\n\nthe molecular geometry must be specified in a separate block:");
208 for (
auto w : what) {
210 else if (
w==
"reference")
reference->param.print(
"dft");
221 return value(
calc->molecule.get_all_coords());
274 return (model==
"ocep") or (model==
"dcep") or (model==
"mrks");
278 return (model==
"dcep");
282 return (model==
"mrks");
287 for (
long i = orbens.
size() - 1; i >= 0; i--) {
288 printf(
" e%2.2lu = %12.8f Eh\n", i, orbens(i));
306 double Drho = rho_diff.
trace();
314 K.set_bra_and_ket(
R_square * nemo, nemo);
332 std::vector<real_function_3d> args={
R_square*rho,numerator,rho,zero,one,lra};
360 for (
int idim=0; idim<3; ++idim) {
365 std::vector<real_function_3d> dnemo=
apply(
world,
D,nemo_copy);
366 u1dnterm+=U1[idim]*
dot(
world,nemo_copy,dnemo);
383 std::vector<vecfuncT> grad_nemo(nemo.size());
384 for (
size_t i = 0; i < nemo.size(); i++) {
389 else grad_nemo[i] =
grad(nemo_copy[i]);
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));
418 if (fockKS.
normf()<1.e-10) {
424 auto [eval, evec] =
syev(fockKS);
425 double homoKS = -eval.max();
427 auto [eval1, evec1] =
syev(fockHF);
428 double homoHF = -eval1.max();
430 double longrange=homoHF-homoKS;
433 for (
int i=0; i<fock1.
dim(0); ++i) fock1(i,i)-=longrange;
446 std::vector<real_function_3d> args={densityKS,ocep_numerator_HF,densityHF,
447 numeratorKS,densityKS,lra_func};
482 std::vector<real_function_3d> args={densityKS,dcep_numerator_HF,densityHF,
483 numeratorKS,densityKS,lra_func};
509 std::vector<real_function_3d> args={densityKS,mrks_numerator_HF,densityHF,
510 numeratorKS,densityKS,lra_func};
515 op.square_denominator=
true;
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];};
577 double Ex = -1.0*
inner(
world, phi, Kphi).sum();
600 const double E_nuc =
calc->molecule.nuclear_repulsion_energy();
601 double energy = E_kin + E_ext + E_J + E_X + E_nuc;
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);
611 return {
energy,E_kin,E_ext,E_J,E_X};
624 for (
long i = 0; i < eigvals.
size(); i++) {
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();
643 double E_1 = -1.0*(E_J + E_K + 2.0*E_Vx) + E_nuc;
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
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
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