78 MADNESS_EXCEPTION(
"implement make_fock operator for your derived NemoBase class",1);
79 return std::shared_ptr<Fock<double,3>>();
83 std::shared_ptr<NuclearCorrelationFactor>
get_ncf_ptr()
const {
88 template<
typename T, std::
size_t NDIM>
92 if (nemo.size()==0)
return;
95 std::vector<Function<T,NDIM> > mos = (metric.is_initialized()) ? metric*nemo : nemo;
99 std::vector<double> invnorm(norms.size());
100 for (std::size_t i = 0; i < norms.size(); ++i)
101 invnorm[i] = 1.0 / norms[i];
108 for (
int i=0; i<s.dim(0); ++i)
Q(i,i) += 1.5;
113 template<
typename T, std::
size_t NDIM>
118 if (nemo.size()==0)
return;
122 std::vector<Function<T,NDIM> > Rnemo = (metric.is_initialized()) ? metric*nemo : nemo;
125 for (
int i=0; i<
Q.dim(0); ++i)
126 for (
int j=0; j<i; ++j)
138 template<
typename T, std::
size_t NDIM>
146 if (not
ncf) need=
true;
158 const std::shared_ptr<PotentialManager> pm,
159 const std::pair<std::string,double> ncf_parameter) {
182 template<
typename T, std::
size_t NDIM>
195 double ke1=
inner(dens,U1dotU1);
204 const std::vector<Function<T,NDIM> > dnemo =
apply(
world,
D, nemo);
222 double ke=2.0*(ke1+ke2+ke3);
230 template<
typename T, std::
size_t NDIM>
234 for (
int i=0; i<nemo.size(); ++i) {
236 ke+=2.0*fnorm2*fnorm2;
238 timer1.
end(
"compute_kinetic_energy1");
247 template<
typename T, std::
size_t NDIM>
251 for (
int i=0; i<
NDIM; ++i) {
252 std::vector< std::shared_ptr< Derivative<T,NDIM> > >
grad=
253 gradient_operator<T,NDIM>(
world);
255 ke+=2.0*fnorm2*fnorm2;
257 timer1.
end(
"compute_kinetic_energy1a");
265 template<
typename T, std::
size_t NDIM>
282 const std::vector<double> oldenergies,
const double bsh_norm,
284 const double econv,
const double dconv)
const {
286 double maxenergychange=fabs(energies.size()-oldenergies.size());
287 for (
auto iter1=energies.begin(), iter2=oldenergies.begin();
288 (iter1!=energies.end() and iter2!=oldenergies.end()); iter1++, iter2++) {
289 maxenergychange=
std::max(maxenergychange,fabs(*iter1 - *iter2));
291 double delta_energy=fabs(energies[0]-oldenergies[0]);
293 bool bsh_conv=
param.converge_bsh_residual() ? bsh_norm<dconv :
true;
294 bool total_energy_conv=
param.converge_total_energy() ? delta_energy<econv :
true;
295 bool each_energy_conv=
param.converge_each_energy() ? maxenergychange<econv*3.0 :
true;
296 bool density_conv=
param.converge_density() ? delta_density<dconv :
true;
299 std::stringstream line;
300 line <<
"convergence: bshresidual, energy change, max energy change, density change "
301 << std::scientific << std::setprecision(1)
302 << bsh_norm <<
" " << delta_energy <<
" "
303 << maxenergychange <<
" " << delta_density;
307 return (bsh_conv and density_conv and each_energy_conv and total_energy_conv);
313 std::shared_ptr<NuclearCorrelationFactor>
ncf;
348 initialize<std::pair<std::string,double> > (
"ncf",{
"slater",2.0},
"nuclear correlation factor");
349 initialize<bool> (
"hessian",
false,
"compute the hessian matrix");
350 initialize<bool> (
"read_cphf",
false,
"read the converged orbital response for nuclear displacements from file");
351 initialize<bool> (
"restart_cphf",
false,
"read the guess orbital response for nuclear displacements from file");
352 initialize<bool> (
"purify_hessian",
false,
"symmetrize the hessian matrix based on atomic charges");
356 std::pair<std::string,double>
ncf()
const {
return get<std::pair<std::string,double> >(
"ncf");}
357 bool hessian()
const {
return get<bool>(
"hessian");}
372 std::string
name()
const {
return "nemo";}
377 print(
"The nemo code computes Hartree-Fock and DFT energies, gradients and hessians using a nuclear correlation factor");
378 print(
"that regularizes the singular nuclear potential. SCF orbitals for the basis for post-SCF calculations like");
379 print(
"excitation energies (cis), correlation energies (cc2), local potentials (oep), etc\n");
380 print(
"A nemo calculation input is mostly identical to a moldft calculation input, but it uses the additional input");
381 print(
"parameter ncf (nuclear correlation factor)\n");
382 print(
"You can print all available calculation parameters by running\n");
383 print(
"nemo --print_parameters\n");
384 print(
"You can perform a simple calculation by running\n");
385 print(
"nemo --geometry=h2o.xyz\n");
386 print(
"provided you have an xyz file in your directory.");
392 print(
"default parameters for the nemo program are");
394 print(
"\n\nthe molecular geometry must be specified in a separate block:");
444 const SCFProtocol&
p,
const std::string& xc_data)
const;
496 return calc->molecule;
517 const int axis)
const;
560 int k=
f.get_impl()->get_k();
625 timer1.
end(
"reproject ncf");
629 poisson = std::shared_ptr<real_convolution_3d>(
674 const std::vector<vecfuncT>&
xi)
const;
701 return calc->vtol / std::min(30.0,
double(
get_calc()->amo.size()));
704 template<
typename solverT>
706 int lo,
int nfunc)
const;
712 coord_3d start(0.0); start[0]=-width;
718 template<
typename T,
size_t NDIM>
722 template<
typename T,
size_t NDIM>
728 template<
typename solverT>
730 int lo,
int nfunc)
const {
731 std::vector < std::vector<Function<double, 3> > > &ulist = solver.get_ulist();
732 std::vector < std::vector<Function<double, 3> > > &rlist = solver.get_rlist();
733 for (
unsigned int iter = 0; iter < ulist.size(); ++iter) {
742 for (
int i=0; i<nfunc; i++) {
751 template<
typename T,
size_t NDIM>
760 template<
typename T,
size_t NDIM>
767 for (std::size_t i=0; i<fsize; ++i) ar &
f[i];
double guess(const coord_3d &r)
Definition: 3dharmonic.cc:127
double w(double t, double eps)
Definition: DKops.h:22
solution protocol for SCF calculations
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition: derivative.h:266
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
static double get_cell_min_width()
Returns the minimum width of any user cell dimension.
Definition: funcdefaults.h:478
A multiresolution adaptive numerical function.
Definition: mra.h:122
double thresh() const
Returns value of truncation threshold. No communication.
Definition: mra.h:567
void set_thresh(double value, bool fence=true)
Sets the value of the truncation threshold. Optional global fence.
Definition: mra.h:577
void clear(bool fence=true)
Clears the function as if constructed uninitialized. Optional fence.
Definition: mra.h:847
bool is_initialized() const
Returns true if the function is initialized.
Definition: mra.h:147
Definition: molecule.h:124
static void print_parameters()
Definition: molecule.cc:110
std::shared_ptr< NuclearCorrelationFactor > ncf
the nuclear correlation factor
Definition: nemo.h:313
virtual std::shared_ptr< Fock< double, 3 > > make_fock_operator() const
Definition: nemo.h:77
double compute_kinetic_energy1a(const std::vector< Function< T, NDIM > > &nemo) const
compute kinetic energy as square of the "analytical" derivative of the orbitals
Definition: nemo.h:248
double compute_kinetic_energy1(const std::vector< Function< T, NDIM > > &nemo) const
compute kinetic energy as square of the "analytical" derivative of the orbitals
Definition: nemo.h:231
std::shared_ptr< NuclearCorrelationFactor > get_ncf_ptr() const
create an instance of the derived object based on the input parameters
Definition: nemo.h:83
virtual void invalidate_factors_and_potentials()
Definition: nemo.h:151
static void normalize(std::vector< Function< T, NDIM > > &nemo, const Function< double, NDIM > metric=Function< double, NDIM >())
normalize the nemos
Definition: nemo.h:89
NemoBase(World &w)
Definition: nemo.h:73
static Tensor< T > Q2(const Tensor< T > &s)
Definition: nemo.h:106
real_function_3d R
the nuclear correlation factor
Definition: nemo.h:316
void construct_nuclear_correlation_factor(const Molecule &molecule, const std::shared_ptr< PotentialManager > pm, const std::pair< std::string, double > ncf_parameter)
Definition: nemo.h:157
virtual bool need_recompute_factors_and_potentials(const double thresh) const
Definition: nemo.h:143
real_function_3d R_square
the square of the nuclear correlation factor
Definition: nemo.h:319
Tensor< double > compute_gradient(const real_function_3d &rhonemo, const Molecule &molecule) const
compute the nuclear gradients
Definition: madness/chem/nemo.cc:94
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
double compute_kinetic_energy2(const std::vector< Function< T, NDIM > > &nemo) const
compute kinetic energy as direct derivative of the orbitals (probably imprecise)
Definition: nemo.h:266
World & world
Definition: nemo.h:310
void orthonormalize(std::vector< Function< T, NDIM > > &nemo, const Function< double, NDIM > metric=Function< double, NDIM >(), const double trantol=FunctionDefaults< NDIM >::get_thresh() *0.01) const
orthonormalize the vectors
Definition: nemo.h:114
virtual ~NemoBase()
Definition: nemo.h:75
bool check_convergence(const std::vector< double > energies, const std::vector< double > oldenergies, const double bsh_norm, const double delta_density, const CalculationParameters ¶m, const double econv, const double dconv) const
Definition: nemo.h:281
Function< typename Tensor< T >::scalar_type, NDIM > compute_density(const std::vector< Function< T, NDIM > > nemo) const
Definition: nemo.h:139
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
bool do_symmetry() const
Definition: nemo.h:692
virtual double value()
Definition: nemo.h:398
projector_irrep get_symmetry_projector() const
return the symmetry_projector
Definition: nemo.h:579
tensorT compute_fock_matrix(const vecfuncT &nemo, const tensorT &occ) const
compute the Fock matrix from scratch
Definition: madness/chem/nemo.cc:291
double trantol() const
return the threshold for vanishing elements in orbital rotations
Definition: nemo.h:700
Nemo(World &world, const commandlineparser &parser)
ctor
Definition: madness/chem/nemo.cc:133
Molecule & molecule() const
return a reference to the molecule
Definition: nemo.h:495
vecfuncT compute_cphf_parallel_term(const size_t iatom, const int iaxis) const
Definition: madness/chem/nemo.cc:1462
void load_function(std::vector< Function< T, NDIM > > &f, const std::string name) const
load a function
Definition: nemo.h:761
void rotate_subspace(World &world, const tensorT &U, solverT &solver, int lo, int nfunc) const
rotate the KAIN subspace (cf. SCF.cc)
Definition: nemo.h:729
Tensor< double > purify_hessian(const Tensor< double > &hessian) const
purify and symmetrize the hessian
Definition: madness/chem/nemo.cc:966
std::shared_ptr< real_convolution_3d > poperatorT
Definition: nemo.h:327
real_function_3d get_coulomb_potential(const vecfuncT &psi) const
return the Coulomb potential
Definition: madness/chem/nemo.cc:589
std::vector< double > compute_energy_regularized(const vecfuncT &nemo, const vecfuncT &Jnemo, const vecfuncT &Knemo, const vecfuncT &Unemo) const
given nemos, compute the HF energy using the regularized expressions for T and V
Definition: madness/chem/nemo.cc:443
Tensor< double > gradient(const Tensor< double > &x)
compute the nuclear gradients
Definition: madness/chem/nemo.cc:780
Tensor< double > hessian(const Tensor< double > &x)
returns the molecular hessian matrix at structure x
Definition: madness/chem/nemo.cc:809
double solve(const SCFProtocol &proto)
solve the HF equations
Definition: madness/chem/nemo.cc:314
Tensor< double > make_incomplete_hessian() const
compute the incomplete hessian
Definition: madness/chem/nemo.cc:1008
std::vector< vecfuncT > compute_all_cphf()
solve the CPHF equation for all displacements
Definition: madness/chem/nemo.cc:1295
Tensor< double > make_incomplete_hessian_response_part(const std::vector< vecfuncT > &xi) const
compute the complementary incomplete hessian
Definition: madness/chem/nemo.cc:1062
real_function_3d make_density(const Tensor< double > &occ, const vecfuncT &nemo) const
make the density (alpha or beta)
Definition: madness/chem/nemo.cc:595
NemoCalculationParameters param
Definition: nemo.h:571
Molecule & molecule()
return a reference to the molecule
Definition: nemo.h:492
std::shared_ptr< SCF > get_calc() const
Definition: nemo.h:482
real_function_3d make_laplacian_density(const real_function_3d &rhonemo) const
the Laplacian of the density
Definition: madness/chem/nemo.cc:637
static void smoothen(real_function_3d &f)
smooth a function by projecting it onto k-1 and then average with k
Definition: nemo.h:559
PCM get_pcm() const
Definition: nemo.h:486
Tensor< double > compute_IR_intensities(const Tensor< double > &normalmodes, const vecfuncT &dens_pt) const
compute the IR intensities in the double harmonic approximation
Definition: madness/chem/nemo.cc:1480
virtual std::shared_ptr< Fock< double, 3 > > make_fock_operator() const
construct the fock operator based on the calculation parameters (K or XC?)
Definition: madness/chem/nemo.cc:247
std::shared_ptr< SCF > calc
Definition: nemo.h:568
static void print_parameters()
Definition: nemo.h:390
bool provides_gradient() const
Override this to return true if the derivative is implemented.
Definition: nemo.h:405
real_function_3d kinetic_energy_potential(const vecfuncT &nemo) const
Definition: madness/chem/nemo.cc:680
void compute_nemo_potentials(const vecfuncT &nemo, vecfuncT &Jnemo, vecfuncT &Knemo, vecfuncT &xcnemo, vecfuncT &pcmnemo, vecfuncT &Unemo) const
compute the reconstructed orbitals, and all potentials applied on nemo
Definition: madness/chem/nemo.cc:516
const NemoCalculationParameters & get_param() const
Definition: nemo.h:484
PCM pcm
polarizable continuum model
Definition: nemo.h:612
void save_function(const std::vector< Function< T, NDIM > > &f, const std::string name) const
save a function
Definition: nemo.h:752
void make_plots(const real_function_3d &f, const std::string &name="function") const
Definition: nemo.h:709
AC< 3 > get_ac() const
Definition: nemo.h:690
bool do_ac() const
Definition: nemo.h:688
projector_irrep symmetry_projector
Definition: nemo.h:574
vecfuncT solve_cphf(const size_t iatom, const int iaxis, const Tensor< double > fock, const vecfuncT &guess, const vecfuncT &rhsconst, const Tensor< double > incomplete_hessian, const vecfuncT ¶llel, const SCFProtocol &p, const std::string &xc_data) const
solve the CPHF equations for the nuclear displacements
Definition: madness/chem/nemo.cc:1141
double coords_sum
sum of square of coords at last solved geometry
Definition: nemo.h:586
std::shared_ptr< real_convolution_3d > poisson
a poisson solver
Definition: nemo.h:590
real_function_3d make_sigma(const real_function_3d &rho1, const real_function_3d &rho2) const
compute the reduced densities sigma (gamma) for GGA functionals
Definition: madness/chem/nemo.cc:751
real_function_3d make_ddensity(const real_function_3d &rhonemo, const int axis) const
make the derivative of the density
Definition: madness/chem/nemo.cc:621
AC< 3 > ac
asymptotic correction for DFT
Definition: nemo.h:593
vecfuncT localize(const vecfuncT &nemo, const double dconv, const bool randomize) const
localize the nemo orbitals
Definition: madness/chem/nemo.cc:233
bool do_pcm() const
Definition: nemo.h:686
bool selftest()
Definition: nemo.h:373
vecfuncT make_cphf_constant_term(const size_t iatom, const int iaxis, const vecfuncT &R2nemo, const real_function_3d &rhonemo) const
compute the constant term for the CPHF equations
Definition: madness/chem/nemo.cc:1090
std::string name() const
Definition: nemo.h:372
static void help()
Definition: nemo.h:375
bool is_dft() const
Definition: nemo.h:684
functor for a local U1 dot U1 potential
Definition: correlationfactor.h:532
interface class to the PCMSolver library
Definition: pcm.h:52
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
class implementing properties of QC models
Definition: QCPropertyInterface.h:11
struct for running a protocol of subsequently tightening precision
Definition: SCFProtocol.h:47
A tensor is a multidimension array.
Definition: tensor.h:317
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
An archive for storing local or parallel data wrapping a BinaryFstreamOutputArchive.
Definition: parallel_archive.h:321
Definition: pointgroupsymmetry.h:98
std::string get_pointgroup() const
get the point group name
Definition: pointgroupsymmetry.h:170
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
Defines/implements plotting interface for functions.
double psi(const Vector< double, 3 > &r)
Definition: hatom_energy.cc:78
static const double v
Definition: hatom_sf_dirac.cc:20
Implements (2nd generation) static load/data balancing for functions.
#define max(a, b)
Definition: lda.h:51
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
optimize the geometrical structure of a molecule
Main include file for MADNESS and defines Function interface.
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
std::shared_ptr< NuclearCorrelationFactor > create_nuclear_correlation_factor(World &world, const Molecule &molecule, const std::shared_ptr< PotentialManager > potentialmanager, const std::string inputline)
create and return a new nuclear correlation factor
Definition: correlationfactor.cc:45
void print_header2(const std::string &s)
medium section heading
Definition: print.cc:54
response_space scale(response_space a, double b)
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
Function< double, NDIM > abssq(const Function< double_complex, NDIM > &z, bool fence=true)
Returns a new function that is the square of the absolute value of the input.
Definition: mra.h:2713
void truncate(World &world, response_space &v, double tol, bool fence)
Definition: basic_operators.cc:30
Function< T, NDIM > project(const Function< T, NDIM > &other, int k=FunctionDefaults< NDIM >::get_k(), double thresh=FunctionDefaults< NDIM >::get_thresh(), bool fence=true)
Definition: mra.h:2399
void set_thresh(World &world, std::vector< Function< T, NDIM > > &v, double thresh, bool fence=true)
Sets the threshold in a vector of functions.
Definition: vmra.h:1251
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:851
void plot_plane(World &world, const Function< double, NDIM > &function, const std::string name)
Definition: funcplot.h:621
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
NDIM & f
Definition: mra.h:2416
std::vector< Function< T, NDIM > > grad(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
shorthand gradient operator
Definition: vmra.h:1818
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
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
vector< functionT > vecfuncT
Definition: corepotential.cc:58
void plot_line(World &world, const char *filename, int npt, const Vector< double, NDIM > &lo, const Vector< double, NDIM > &hi, const opT &op)
Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi.
Definition: funcplot.h:438
std::string name(const FuncType &type, const int ex=-1)
Definition: ccpairfunction.h:28
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 matrix_inner(DistributedMatrix< T > &A, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition: distpm.cc:46
std::vector< double > norm2s(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norms of a vector of functions.
Definition: vmra.h:827
static long abs(long a)
Definition: tensor.h:218
Implementation of Krylov-subspace nonlinear equation solver.
Implements most functionality of separated operators.
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
Definition: test_ar.cc:204
Definition: CalculationParameters.h:51
std::string ac_data() const
Definition: CalculationParameters.h:196
int print_level() const
Definition: CalculationParameters.h:188
double lo() const
Definition: CalculationParameters.h:177
std::string pcm_data() const
Definition: CalculationParameters.h:195
Definition: molecular_optimizer.h:46
virtual Molecule & molecule()
return the molecule of the target
Definition: molecular_optimizer.h:49
class holding parameters for a nemo calculation beyond the standard dft parameters from moldft
Definition: nemo.h:333
void initialize_nemo_parameters()
Definition: nemo.h:347
std::pair< std::string, double > ncf() const
Definition: nemo.h:356
bool hessian() const
Definition: nemo.h:357
NemoCalculationParameters(World &world, const commandlineparser &parser)
Definition: nemo.h:335
NemoCalculationParameters(const CalculationParameters ¶m)
Definition: nemo.h:339
NemoCalculationParameters()
Definition: nemo.h:343
very simple command line parser
Definition: commandlineparser.h:15
Definition: timing_utilities.h:9
double end(const std::string msg)
Definition: timing_utilities.h:56
Definition: dirac-hatom.cc:108
InputParameters param
Definition: tdse.cc:203
static const std::size_t NDIM
Definition: testpdiff.cc:42
std::size_t axis
Definition: testpdiff.cc:59
Defines operations on vectors of Functions.