40#ifndef MADNESS_CHEM_SCF_H__INCLUDED
41#define MADNESS_CHEM_SCF_H__INCLUDED
62#include <madness/tensor/tensor_json.hpp>
67typedef std::shared_ptr<WorldDCPmapInterface<Key<3> > >
pmapT;
69typedef std::shared_ptr<FunctionFunctorInterface<double, 3> >
functorT;
71typedef std::vector<functionT>
vecfuncT;
84template<
typename T,
int NDIM>
92 if (key.
level() < 1) {
113 x = (x * x * (3. - 2. * x));
120 double x = rsim[0], y = rsim[1], z = rsim[2];
121 double lo = 0.0625, hi = 1.0 -
lo, result = 1.0;
122 double rlo = 1.0 /
lo;
125 result *=
mask1(x * rlo);
127 result *=
mask1((1.0 - x) * rlo);
129 result *=
mask1(y * rlo);
131 result *=
mask1((1.0 - y) * rlo);
133 result *=
mask1(z * rlo);
135 result *=
mask1((1.0 - z) * rlo);
141class DipoleFunctor :
public FunctionFunctorInterface<double, 3> {
163 double xi = 1.0, yj = 1.0, zk = 1.0;
164 for (
int p = 0;
p <
i; ++
p)
xi *= r[0];
165 for (
int p = 0;
p <
j; ++
p) yj *= r[1];
166 for (
int p = 0;
p <
k; ++
p) zk *= r[2];
173 std::map<std::string, std::vector<double>>
e_data;
183 void add_data(std::map<std::string, double> values);
218 std::vector<std::shared_ptr<real_derivative_3d> >
gradop;
236 print(
"The moldft code computes Hartree-Fock and DFT energies and gradients, It is the fastest code in MADNESS");
237 print(
"and considered the reference implementation. No nuclear correlation factor can be used");
238 print(
"SCF orbitals are the basis for post-SCF calculations like");
239 print(
"excitation energies (cis), correlation energies (cc2), local potentials (oep), etc\n\n");
240 print(
"You can print all available calculation parameters by running\n");
241 print(
"moldft --print_parameters\n");
242 print(
"You can perform a simple calculation by running\n");
243 print(
"moldft --geometry=h2o.xyz\n");
244 print(
"provided you have an xyz file in your directory.");
250 print(
"default parameters for the moldft program are");
252 print(
"\n\nthe molecular geometry must be specified in a separate block:");
258 template<std::
size_t NDIM>
264 else if (
thresh >= 0.9e-4)
266 else if (
thresh >= 0.9e-6)
268 else if (
thresh >= 0.9e-8)
294 gradop = gradient_operator<double, 3>(world);
298 for (
int i = 0; i < 3; ++i) (*
gradop[i]).set_bspline1();
300 for (
int i = 0; i < 3; ++i) (*
gradop[i]).set_ble1();
353 const tensorT& occ,
const int nmo)
const;
389 const functionT& vlocal,
double& exc,
double& enl,
int ispin);
400 double& maxabsval)
const;
407 double& ekinetic)
const;
431 const double thresh_degenerate)
const;
455 int lo,
int nfunc,
double trantol)
const;
458 int lo,
int nfunc,
double trantol)
const;
464 double& bsh_residual,
double& update_residual);
476 vecfuncT& mo_new, std::string spin)
const;
529 std::string
name()
const {
return "Molecularenerg"; }
536 double xsq = x.
sumsq();
556 bool have_initial_guess =
false;
560 have_initial_guess =
true;
567 if (not have_initial_guess) {
586 int nvalpha_start, nv_old;
589 if (proto == 0 && nvalpha > 0) {
592 nvalpha_start = nvalpha;
595 nv_old = nvalpha_start;
597 for (
int nv = nvalpha_start; nv >= nvalpha; nv -= nvalpha) {
599 if (nv > 0 &&
world.
rank() == 0) std::cout <<
"Running with " << nv <<
" virtual states" << std::endl;
604 if (nvbeta == nvalpha) {
666 rho.
gaxpy(1.0, brho, 1.0);
679 rho.
gaxpy(1.0, brho, 1.0);
686 nlohmann::json j = {};
687 vec_pair_ints int_vals;
688 vec_pair_T<double> double_vals;
689 vec_pair_tensor_T<double> double_tensor_vals;
693 nlohmann::json calc_precision={ };
702 int_vals.push_back({
"calcinfo_nmo",
param.nmo_alpha() +
param.nmo_beta()});
703 int_vals.push_back({
"calcinfo_nalpha",
param.nalpha()});
704 int_vals.push_back({
"calcinfo_nbeta",
param.nbeta()});
711 double_tensor_vals.push_back({
"scf_eigenvalues_a",
calc.
aeps});
712 if (
param.nbeta() != 0 && !
param.spin_restricted()) {
713 double_tensor_vals.push_back({
"scf_eigenvalues_b",
calc.
beps});
716 to_json(j, double_tensor_vals);
720 j[
"precision"]=calc_precision;
721 j[
"molecule"]=mol_json;
Operators for the molecular HF and DFT code.
A MADNESS functor to compute either x, y, or z.
Definition preal.cc:136
Contracted Gaussian basis.
Definition madness/chem/molecularbasis.h:465
void read_file(std::string filename)
read the atomic basis set from file
Definition molecularbasis.cc:118
Provides the common functionality/interface of all 1D convolutions.
Definition convolution1d.h:257
DipoleFunctor(int axis)
Definition SCF.h:145
double operator()(const coordT &x) const
Definition SCF.h:147
const int axis
Definition solver.h:167
Manages data associated with a row/column/block distributed array.
Definition distributed_matrix.h:388
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:204
static int get_k()
Returns the default wavelet order.
Definition funcdefaults.h:266
static void set_apply_randomize(bool value)
Sets the random load balancing for integral operators flag.
Definition funcdefaults.h:397
static void set_thresh(double value)
Sets the default threshold.
Definition funcdefaults.h:286
static void set_k(int value)
Sets the default wavelet order.
Definition funcdefaults.h:273
static const double & get_thresh()
Returns the default threshold.
Definition funcdefaults.h:279
static void set_autorefine(bool value)
Sets the default adaptive autorefinement flag.
Definition funcdefaults.h:363
static void set_project_randomize(bool value)
Sets the random load balancing for projection flag.
Definition funcdefaults.h:408
static void set_initial_level(int value)
Sets the default initial projection level.
Definition funcdefaults.h:303
static void set_cubic_cell(double lo, double hi)
Sets the user cell to be cubic with each dimension having range [lo,hi].
Definition funcdefaults.h:461
static void set_refine(bool value)
Sets the default adaptive refinement flag.
Definition funcdefaults.h:351
FunctionFactory implements the named-parameter idiom for Function.
Definition function_factory.h:86
Abstract base class interface required for functors used as input to Functions.
Definition function_interface.h:68
FunctionNode holds the coefficients, etc., at each node of the 2^NDIM-tree.
Definition funcimpl.h:124
bool is_leaf() const
Returns true if this does not have children.
Definition funcimpl.h:210
A multiresolution adaptive numerical function.
Definition mra.h:122
void broaden(const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), bool fence=true) const
Inplace broadens support in scaling function basis.
Definition mra.h:836
void sum_down(bool fence=true) const
Sums scaling coeffs down tree restoring state with coeffs only at leaves. Optional fence....
Definition mra.h:798
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition mra.h:1099
const Function< T, NDIM > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition mra.h:775
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition mra.h:981
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:66
Level level() const
Definition key.h:159
bool selftest()
Definition SCF.h:531
bool provides_gradient() const
Override this to return true if the derivative is implemented.
Definition SCF.h:533
double value(const Tensor< double > &x)
Should return the value of the objective function.
Definition SCF.h:535
World & world
Definition SCF.h:521
SCF & calc
Definition SCF.h:522
double coords_sum
Definition SCF.h:523
void output_calc_info_schema()
Definition SCF.h:685
madness::Tensor< double > gradient(const Tensor< double > &x)
Should return the derivative of the function.
Definition SCF.h:659
void energy_and_gradient(const Molecule &molecule, double &energy, Tensor< double > &gradient)
Definition SCF.h:672
std::string name() const
Definition SCF.h:529
MolecularEnergy(World &world, SCF &calc)
Definition SCF.h:526
Definition molecule.h:124
void set_all_coords(const madness::Tensor< double > &newcoords)
Definition molecule.cc:419
madness::Tensor< double > get_all_coords() const
Definition molecule.cc:397
size_t natom() const
Definition molecule.h:387
static void print_parameters()
Definition molecule.cc:110
json to_json() const
Definition molecule.cc:457
GeometryParameters parameters
Definition molecule.h:257
A MADNESS functor to compute the cartesian moment x^i * y^j * z^k (i, j, k integer and >= 0)
Definition SCF.h:154
const int j
Definition SCF.h:156
double operator()(const coordT &r) const
Definition SCF.h:162
MomentFunctor(int i, int j, int k)
Definition SCF.h:158
const int i
Definition SCF.h:156
MomentFunctor(const std::vector< int > &x)
Definition SCF.h:160
const int k
Definition SCF.h:156
interface class to the PCMSolver library
Definition pcm.h:52
void set_user_defined_value(const std::string &key, const T &value)
Definition QCCalculationParametersBase.h:533
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
class implementing properties of QC models
Definition QCPropertyInterface.h:11
void copy_data(World &world, const SCF &other)
Definition SCF.cc:304
void do_plots(World &world)
Definition SCF.cc:487
tensorT derivatives(World &world, const functionT &rho) const
Definition SCF.cc:1389
void save_mos(World &world)
Definition SCF.cc:318
static vecfuncT project_ao_basis_only(World &world, const AtomicBasisSet &aobasis, const Molecule &molecule)
Definition SCF.cc:581
const tensorT & get_aocc() const
getter for the occupation numbers, alpha spin
Definition SCF.h:317
std::shared_ptr< GTHPseudopotential< double > > gthpseudopotential
Definition SCF.h:190
void make_nuclear_potential(World &world)
Definition SCF.cc:556
std::vector< int > aset
Definition SCF.h:205
static void print_parameters()
Definition SCF.h:248
AtomicBasisSet aobasis
Definition SCF.h:195
vecfuncT ao
MRA projection of the minimal basis set.
Definition SCF.h:208
scf_data e_data
Definition SCF.h:198
void initial_guess(World &world)
Definition SCF.cc:735
vecfuncT apply_potential(World &world, const tensorT &occ, const vecfuncT &amo, const functionT &vlocal, double &exc, double &enl, int ispin)
Definition SCF.cc:1309
vecfuncT amo
alpha and beta molecular orbitals
Definition SCF.h:201
vecfuncT compute_residual(World &world, tensorT &occ, tensorT &fock, const vecfuncT &psi, vecfuncT &Vpsi, double &err)
Definition SCF.cc:1500
const tensorT & get_bocc() const
getter for the occupation numbers, alpha spin
Definition SCF.h:320
std::vector< int > at_to_bf
Definition SCF.h:210
std::vector< int > bset
Definition SCF.h:205
distmatT kinetic_energy_matrix(World &world, const vecfuncT &v) const
Definition SCF.cc:653
tensorT bocc
Definition SCF.h:213
double vtol
Definition SCF.h:219
tensorT dipole(World &world, const functionT &rho) const
compute the total dipole moment of the molecule
Definition SCF.cc:1462
PCM pcm
Definition SCF.h:194
bool is_spin_restricted() const
Definition SCF.h:322
poperatorT coulop
Definition SCF.h:217
double make_dft_energy(World &world, const vecfuncT &vf, int ispin)
Definition SCF.h:382
void output_calc_info_schema() const
Definition SCF.cc:157
static void help()
Definition SCF.h:234
void update_subspace(World &world, vecfuncT &Vpsia, vecfuncT &Vpsib, tensorT &focka, tensorT &fockb, subspaceT &subspace, tensorT &Q, double &bsh_residual, double &update_residual)
Definition SCF.cc:1776
Molecule molecule
Definition SCF.h:191
void load_mos(World &world)
Definition SCF.cc:358
tensorT beps
Definition SCF.h:216
std::vector< poperatorT > make_bsh_operators(World &world, const tensorT &evals) const
Definition SCF.cc:1278
tensorT diag_fock_matrix(World &world, tensorT &fock, vecfuncT &psi, vecfuncT &Vpsi, tensorT &evals, const tensorT &occ, const double thresh) const
diagonalize the fock matrix, taking care of degenerate states
Definition SCF.cc:1674
const vecfuncT & get_amo() const
getter for the molecular orbitals, alpha spin
Definition SCF.h:311
static void analyze_vectors(World &world, const vecfuncT &mo, const vecfuncT &ao, double vtol, const Molecule &molecule, const int print_level, const AtomicBasisSet &aobasis, const tensorT &occ=tensorT(), const tensorT &energy=tensorT(), const std::vector< int > &set=std::vector< int >())
Definition SCF.cc:595
XCfunctional xc
Definition SCF.h:193
void initial_load_bal(World &world)
Definition SCF.cc:1220
std::vector< int > at_nbf
Definition SCF.h:210
double do_step_restriction(World &world, const vecfuncT &mo, vecfuncT &mo_new, std::string spin) const
perform step restriction following the KAIN solver
Definition SCF.cc:1918
std::vector< std::shared_ptr< real_derivative_3d > > gradop
Definition SCF.h:218
void set_print_timings(const bool value)
Definition SCF.cc:300
tensorT get_fock_transformation(World &world, const tensorT &overlap, tensorT &fock, tensorT &evals, const tensorT &occ, const double thresh_degenerate) const
compute the unitary transformation that diagonalizes the fock matrix
Definition SCF.cc:1643
bool restart_aos(World &world)
Definition SCF.cc:683
tensorT aeps
orbital energies for alpha and beta orbitals
Definition SCF.h:216
functionT make_coulomb_potential(const functionT &rho) const
make the Coulomb potential given the total density
Definition SCF.h:410
tensorT make_fock_matrix(World &world, const vecfuncT &psi, const vecfuncT &Vpsi, const tensorT &occ, double &ekinetic) const
Definition SCF.cc:1555
Tensor< double > twoint(World &world, const vecfuncT &psi) const
Compute the two-electron integrals over the provided set of orbitals.
Definition SCF.cc:1613
tensorT aocc
occupation numbers for alpha and beta orbitals
Definition SCF.h:213
void set_protocol(World &world, double thresh)
Definition SCF.h:259
static functionT make_lda_potential(World &world, const functionT &arho)
Definition SCF.cc:1301
functionT make_density(World &world, const tensorT &occ, const vecfuncT &v) const
Definition SCF.cc:1237
void rotate_subspace(World &world, const tensorT &U, subspaceT &subspace, int lo, int nfunc, double trantol) const
Definition SCF.cc:1742
void reset_aobasis(const std::string &aobasisname)
Definition SCF.h:341
void project(World &world)
Definition SCF.cc:534
double converged_for_thresh
mos are converged for this threshold
Definition SCF.h:221
CalculationParameters param
Definition SCF.h:192
double current_energy
Definition SCF.h:220
void loadbal(World &world, functionT &arho, functionT &brho, functionT &arho_old, functionT &brho_old, subspaceT &subspace)
Definition SCF.cc:1707
void vector_stats(const std::vector< double > &v, double &rms, double &maxabsval) const
Definition SCF.cc:1488
vecfuncT project_ao_basis(World &world, const AtomicBasisSet &aobasis)
Definition SCF.cc:573
std::vector< int > group_orbital_sets(World &world, const tensorT &eps, const tensorT &occ, const int nmo) const
group orbitals into sets of similar orbital energies for localization
Definition SCF.cc:1195
complex_functionT APPLY(const complex_operatorT *q1d, const complex_functionT &psi)
Definition SCF.h:486
functionT mask
Definition SCF.h:196
std::shared_ptr< PotentialManager > potentialmanager
Definition SCF.h:189
const vecfuncT & get_bmo() const
getter for the molecular orbitals, beta spin
Definition SCF.h:314
vecfuncT bmo
Definition SCF.h:201
void solve(World &world)
Definition SCF.cc:2040
void orthonormalize(World &world, vecfuncT &amo_new) const
orthonormalize the vectors
Definition SCF.cc:1990
Convolutions in separated form (including Gaussian)
Definition operator.h:136
A tensor is a multidimension array.
Definition tensor.h:317
Tensor< T > reshape(int ndimnew, const long *d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition tensor.h:1384
T sumsq() const
Returns the sum of the squares of the elements.
Definition tensor.h:1669
Tensor< T > flat()
Returns new view/tensor rehshaping to flat (1-d) tensor.
Definition tensor.h:1555
A simple, fixed dimension vector.
Definition vector.h:64
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition worldgop.cc:161
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
WorldGopInterface & gop
Global operations.
Definition world.h:205
Simplified interface to XC functionals.
Definition xcfunctional.h:43
std::map< std::string, std::vector< double > > e_data
Definition SCF.h:173
int iter
Definition SCF.h:174
void to_json(json &j) const
Definition SCF.cc:213
void add_data(std::map< std::string, double > values)
Definition SCF.cc:190
scf_data()
Definition SCF.cc:200
void print_data()
Definition SCF.cc:225
Declaration of core potential related class.
double(* energy)()
Definition derivatives.cc:58
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
static double lo
Definition dirac-hatom.cc:23
double psi(const Vector< double, 3 > &r)
Definition hatom_energy.cc:78
static const double v
Definition hatom_sf_dirac.cc:20
Main include file for MADNESS and defines Function interface.
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
static SeparatedConvolution< double, 3 > * CoulombOperatorPtr(World &world, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition operator.h:1762
void print_header2(const std::string &s)
medium section heading
Definition print.cc:54
Function< TENSOR_RESULT_TYPE(typename opT::opT, R), NDIM > apply_1d_realspace_push(const opT &op, const Function< R, NDIM > &f, int axis, bool fence=true)
Definition mra.h:2231
Vector< double, 3 > coordT
Definition corepotential.cc:54
nlohmann::json json
Definition QCCalculationParametersBase.h:27
static void user_to_sim(const Vector< double, NDIM > &xuser, Vector< double, NDIM > &xsim)
Convert user coords (cell[][]) to simulation coords ([0,1]^ndim)
Definition funcdefaults.h:524
Tensor< double > tensorT
Definition distpm.cc:21
Function< std::complex< double >, 3 > complex_functionT
Definition SCF.h:79
std::vector< pairvecfuncT > subspaceT
Definition SCF.h:73
double mask1(double x)
Definition SCF.h:103
DistributedMatrix< double > distmatT
Definition SCF.h:75
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition SCF.h:72
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
FunctionFactory< double, 3 > factoryT
Definition corepotential.cc:57
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
std::shared_ptr< operatorT> poperatorT
Definition SCF.h:78
Function< double, 3 > functionT
Definition corepotential.cc:56
NDIM & f
Definition mra.h:2416
void to_json(nlohmann::json &j)
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition corepotential.cc:55
std::vector< complex_functionT > cvecfuncT
Definition SCF.h:80
vector< functionT > vecfuncT
Definition corepotential.cc:58
static double mask3(const coordT &ruser)
Definition SCF.h:117
std::shared_ptr< WorldDCPmapInterface< Key< 3 > > > pmapT
Definition SCF.h:67
Convolution1D< double_complex > complex_operatorT
Definition SCF.h:81
SeparatedConvolution< double, 3 > operatorT
Definition SCF.h:77
double Q(double a)
Definition relops.cc:20
static const double thresh
Definition rk.cc:45
static const long k
Definition rk.cc:44
const double xi
Exponent for delta function approx.
Definition siam_example.cc:60
Defines interfaces for optimization and non-linear equation solvers.
Definition CalculationParameters.h:51
std::vector< double > protocol() const
Definition CalculationParameters.h:202
int nv_factor() const
Definition CalculationParameters.h:167
int k() const
Definition CalculationParameters.h:179
int nmo_alpha() const
Definition CalculationParameters.h:169
bool save() const
Definition CalculationParameters.h:203
double dconv() const
Definition CalculationParameters.h:142
double econv() const
Definition CalculationParameters.h:141
int print_level() const
Definition CalculationParameters.h:188
std::string deriv() const
Definition CalculationParameters.h:193
double L() const
Definition CalculationParameters.h:178
std::string aobasis() const
Definition CalculationParameters.h:200
int nalpha() const
Definition CalculationParameters.h:162
int nbeta() const
Definition CalculationParameters.h:163
int nmo_beta() const
Definition CalculationParameters.h:170
bool no_compute() const
Definition CalculationParameters.h:175
bool spin_restricted() const
Definition CalculationParameters.h:174
double lo() const
Definition CalculationParameters.h:177
std::string pcm_data() const
Definition CalculationParameters.h:195
bool restart() const
Definition CalculationParameters.h:204
bool restartao() const
Definition CalculationParameters.h:205
Definition convolution1d.h:849
double eprec() const
Definition molecule.h:228
The interface to be provided by functions to be optimized.
Definition solvers.h:176
very simple command line parser
Definition commandlineparser.h:15
double parent_value
Definition SCF.h:87
lbcost(double leaf_value=1.0, double parent_value=0.0)
Definition SCF.h:89
double leaf_value
Definition SCF.h:86
double operator()(const Key< NDIM > &key, const FunctionNode< T, NDIM > &node) const
Definition SCF.h:91
Class to compute the energy functional.
Definition xcfunctional.h:360
InputParameters param
Definition tdse.cc:203
static const std::size_t NDIM
Definition testpdiff.cc:42
static Molecule molecule
Definition testperiodicdft.cc:38
static Subspace * subspace
Definition testperiodicdft.cc:40