MADNESS  0.10.1
nemo.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30  */
31 
32 /*!
33  \file examples/nemo.h
34  \brief solve the HF equations using numerical exponential MOs
35 
36  The source is
37  <a href=http://code.google.com/p/m-a-d-n-e-s-s/source/browse/local
38  /trunk/src/apps/examples/nemo.h>here</a>.
39 
40  */
41 
42 #ifndef NEMO_H_
43 #define NEMO_H_
44 
45 #include <madness/mra/mra.h>
46 #include <madness/mra/funcplot.h>
47 #include <madness/mra/operator.h>
48 #include <madness/mra/lbdeux.h>
49 #include<madness/chem/SCF.h>
54 #include <madness/mra/nonlinsol.h>
55 #include <madness/mra/vmra.h>
56 #include<madness/chem/pcm.h>
57 #include<madness/chem/AC.h>
62 
63 namespace madness {
64 
65 class PNO;
66 class OEP;
67 
68 
70 
71 public:
72 
74 
75  virtual ~NemoBase() {}
76 
77  virtual std::shared_ptr<Fock<double,3>> make_fock_operator() const {
78  MADNESS_EXCEPTION("implement make_fock operator for your derived NemoBase class",1);
79  return std::shared_ptr<Fock<double,3>>();
80  }
81 
82  /// create an instance of the derived object based on the input parameters
83  std::shared_ptr<NuclearCorrelationFactor> get_ncf_ptr() const {
84  return ncf;
85  }
86 
87  /// normalize the nemos
88  template<typename T, std::size_t NDIM>
89  void static normalize(std::vector<Function<T,NDIM> >& nemo,
91 
92  if (nemo.size()==0) return;
93  World& world=nemo[0].world();
94  // compute the norm of the reconstructed orbitals, includes the factor
95  std::vector<Function<T,NDIM> > mos = (metric.is_initialized()) ? metric*nemo : nemo;
96  std::vector<double> norms = norm2s(world, mos);
97 
98  // scale the nemos, excludes the nuclear correlation factor
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];
102  scale(world, nemo, invnorm);
103  }
104 
105  template<typename T>
106  static Tensor<T> Q2(const Tensor<T>& s) {
107  Tensor<T> Q = -0.5*s;
108  for (int i=0; i<s.dim(0); ++i) Q(i,i) += 1.5;
109  return Q;
110  }
111 
112  /// orthonormalize the vectors
113  template<typename T, std::size_t NDIM>
114  void orthonormalize(std::vector<Function<T,NDIM> >& nemo,
116  const double trantol=FunctionDefaults<NDIM>::get_thresh()*0.01) const {
117 
118  if (nemo.size()==0) return;
119  normalize(nemo,metric);
120  double maxq;
121  do {
122  std::vector<Function<T,NDIM> > Rnemo = (metric.is_initialized()) ? metric*nemo : nemo;
123  Tensor<T> Q = Q2(matrix_inner(world, Rnemo, Rnemo));
124  maxq=0.0;
125  for (int i=0; i<Q.dim(0); ++i)
126  for (int j=0; j<i; ++j)
127  maxq = std::max(maxq,std::abs(Q(i,j)));
128 
129  Q.screen(trantol); // ???? Is this really needed?
130  nemo = transform(world, nemo, Q, trantol, true);
131  truncate(world, nemo);
132 // if (world.rank() == 0) print("ORTHOG2: maxq trantol", maxq, trantol);
133 
134  } while (maxq>0.01);
135  normalize(nemo,metric);
136  }
137 
138  template<typename T, std::size_t NDIM>
140  return sum(world,abssq(world,nemo)).truncate();
141  }
142 
143  virtual bool need_recompute_factors_and_potentials(const double thresh) const {
144  bool need=false;
145  if ((not R.is_initialized()) or (R.thresh()>thresh)) need=true;
146  if (not ncf) need=true;
147  if ((not R_square.is_initialized()) or (R_square.thresh()>thresh)) need=true;
148  return need;
149  };
150 
152  R.clear();
153  R_square.clear();
154  ncf.reset();
155  };
156 
158  const std::shared_ptr<PotentialManager> pm,
159  const std::pair<std::string,double> ncf_parameter) {
160 
161  // construct the nuclear correlation factor:
162  if (not ncf) {
164  }
165 
166  // re-project the ncf
167  ncf->initialize(FunctionDefaults<3>::get_thresh());
168  R = ncf->function();
170  R_square = ncf->square();
172  }
173 
174  /// compute the nuclear gradients
176  const Molecule& molecule) const;
177 
178  /// compute kinetic energy as square of the "analytical" expectation value
179 
180  /// @param[in] the nemo orbitals F
181  /// @return T = 1/2 \sum_i \int R^2 U1.U1 F^2 + 2 R^2 U1.grad(F) + R^2 grad(F)^2
182  template<typename T, std::size_t NDIM>
183  double compute_kinetic_energy(const std::vector<Function<T,NDIM> >& nemo) const {
184 
185  // T = 0.5\sum_i \int R^2 U1.U1 F^2 - 2 R^2 U1.grad(F) F + R^2 grad(F)^2
186  // = 0.5 (<U1.U1 | rho > + <R^2|grad(F)^2> - 2<R^2 | U1.grad(F) >)
187  // note: U1=-grad(R)/R
188  //auto id=nemo.front().world().id();
189  //auto id1=R_square.world().id();
190  //auto worldid=world.id();
191  world.gop.fence();
192  real_function_3d dens=dot(world,nemo,nemo)*R_square;
195  double ke1=inner(dens,U1dotU1);
196 
197  double ke2=0.0;
198  double ke3=0.0;
199  //double ke3_real=0.0;
200  //double ke3_imag=0.0;
201 
202  for (size_t axis = 0; axis < NDIM; axis++) {
203  real_derivative_3d D = free_space_derivative<double, NDIM>(world, axis);
204  const std::vector<Function<T,NDIM> > dnemo = apply(world, D, nemo);
205 
206  real_function_3d term2=dot(world,dnemo,nemo)*ncf->U1(axis);
207  double tmp=-2.0*inner(R_square,term2);
208 // ke3_real -=2.0*std::real(tmp);
209 // ke3_imag -=2.0*std::imag(tmp);
210  ke3 +=tmp;
211 
212  const real_function_3d term1=dot(world,dnemo,dnemo);
213  world.gop.fence();
214  ke2 += inner(term1,R_square);
215 
216  }
217 // if (ke3_imag>1.e-8) {
218 // print("kinetic energy, imaginary part: ",ke3_imag);
219 // MADNESS_EXCEPTION("imaginary kinetic energy",1);
220 // }
221 // double ke=2.0*(ke1+ke2+ke3_real); // closed shell
222  double ke=2.0*(ke1+ke2+ke3); // closed shell
223  return 0.5*ke;
224  }
225 
226  /// compute kinetic energy as square of the "analytical" derivative of the orbitals
227 
228  /// @param[in] the nemo orbitals F
229  /// @return T = 1/2 \sum_i || grad(R)*F_i + R*grad(F_i)||^2
230  template<typename T, std::size_t NDIM>
231  double compute_kinetic_energy1(const std::vector<Function<T,NDIM> >& nemo) const {
232  timer timer1(world);
233  double ke=0.0;
234  for (int i=0; i<nemo.size(); ++i) {
235  double fnorm2=norm2(world,-1.0*R*ncf->U1vec()*nemo[i] + R*grad(nemo[i]));
236  ke+=2.0*fnorm2*fnorm2;
237  }
238  timer1.end("compute_kinetic_energy1");
239  return 0.5*ke;
240  }
241 
242 
243  /// compute kinetic energy as square of the "analytical" derivative of the orbitals
244 
245  /// @param[in] the nemo orbitals F
246  /// @return T = 1/2 \sum_i || grad(R)*F_i + R*grad(F_i)||^2
247  template<typename T, std::size_t NDIM>
248  double compute_kinetic_energy1a(const std::vector<Function<T,NDIM> >& nemo) const {
249  timer timer1(world);
250  double ke=0.0;
251  for (int i=0; i<NDIM; ++i) {
252  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
253  gradient_operator<T,NDIM>(world);
254  double fnorm2=norm2(world,R*(-1.0*ncf->U1(i)*nemo + apply(world,*(grad[i]),nemo)));
255  ke+=2.0*fnorm2*fnorm2;
256  }
257  timer1.end("compute_kinetic_energy1a");
258  return 0.5*ke;
259  }
260 
261  /// compute kinetic energy as direct derivative of the orbitals (probably imprecise)
262 
263  /// @param[in] the nemo orbitals F
264  /// @return T = 1/2 \sum_i || grad(R*F_i)||^2
265  template<typename T, std::size_t NDIM>
266  double compute_kinetic_energy2(const std::vector<Function<T,NDIM> >& nemo) const {
267 
268  // it's ok to use phi here, no regularization necessary for this eigenvalue
269  double E_kin = 0.0;
270  for (int axis = 0; axis < 3; axis++) {
271  real_derivative_3d D = free_space_derivative<double, 3>(world, axis);
272  const vecfuncT dphi = apply(world, D,R*nemo);
273  E_kin += 0.5 * (inner(world, dphi, dphi)).sum();
274  // -1/2 sum <Psi|Nabla^2|Psi> = 1/2 sum <NablaPsi|NablaPsi> (integration by parts)
275  }
276  E_kin *= 2.0; // 2 because closed shell
277  return E_kin;
278  }
279 
280 
281  bool check_convergence(const std::vector<double> energies,
282  const std::vector<double> oldenergies, const double bsh_norm,
283  const double delta_density, const CalculationParameters& param,
284  const double econv, const double dconv) const {
285 
286  double maxenergychange=fabs(energies.size()-oldenergies.size()); // >0 if oldenergyvec not initialized
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));
290  }
291  double delta_energy=fabs(energies[0]-oldenergies[0]);
292 
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;
297 
298  if (world.rank()==0 and param.print_level()>2) {
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;
304  print(line.str());
305  }
306 
307  return (bsh_conv and density_conv and each_energy_conv and total_energy_conv);
308  }
309 
311 
312  /// the nuclear correlation factor
313  std::shared_ptr<NuclearCorrelationFactor> ncf;
314 
315  /// the nuclear correlation factor
317 
318  /// the square of the nuclear correlation factor
320 
321 
322 };
323 
324 
325 /// The Nemo class
326 class Nemo: public NemoBase, public QCPropertyInterface {
327  typedef std::shared_ptr<real_convolution_3d> poperatorT;
328  friend class PNO;
329  friend class TDHF;
330 
331 public:
332  /// class holding parameters for a nemo calculation beyond the standard dft parameters from moldft
334 
337  }
338 
341  }
342 
345  }
346 
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");
353  set_derived_value("k",7);
354  }
355 
356  std::pair<std::string,double> ncf() const {return get<std::pair<std::string,double> >("ncf");}
357  bool hessian() const {return get<bool>("hessian");}
358 
359  };
360 
361 
362 public:
363 
364  /// ctor
365 
366  /// @param[in] world1 the world
367  /// @param[in] calc the SCF
368 // Nemo(World& world1, std::shared_ptr<SCF> calc, const std::string inputfile);
369 
370  Nemo(World& world, const commandlineparser& parser);
371 
372  std::string name() const {return "nemo";}
373  bool selftest() {return false;}
374 
375  static void help() {
376  print_header2("help page for 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.");
387 
388  }
389 
390  static void print_parameters() {
392  print("default parameters for the nemo program are");
393  param.print("dft","end");
394  print("\n\nthe molecular geometry must be specified in a separate block:");
396  }
397 
398  virtual double value() {return value(calc->molecule.get_all_coords());}
399 
400  virtual double value(const Tensor<double>& x);
401 
402  /// compute the nuclear gradients
404 
405  bool provides_gradient() const {return true;}
406 
407  /// returns the molecular hessian matrix at structure x
409 
410  /// construct the fock operator based on the calculation parameters (K or XC?)
411  virtual std::shared_ptr<Fock<double,3>> make_fock_operator() const;
412 
413  /// purify and symmetrize the hessian
414 
415  /// The hessian should be symmetric, but it is not, because
416  /// \f[
417  /// \langle i^{Y_B}|H^{X_A}|i\rangle \neq \langle i|H^{X_A}|i^{Y_B}\rangle
418  /// \f]
419  /// does holds analytically, but not numerically. If the two numbers
420  /// differ, pick the more trustworthy, which is the one with a heavy
421  /// atom causing the perturbed density and the light atom being the
422  /// nuclear singularity.
423  /// @param[in] hessian the raw hessian
424  /// @return a symmetrized hessian
426 
427 
428  /// solve the CPHF equations for the nuclear displacements
429 
430  /// this function computes that part of the orbital response that is
431  /// orthogonal to the occupied space. If no NCF's are used this
432  /// corresponds to the normal response. If NCF's are used the part
433  /// parallel to the occupied space must be added!
434  /// \f[
435  /// F^X = F^\perp + F^\parallel
436  /// \f]
437  /// cf parallel_CPHF()
438  /// @param[in] iatom the atom A to be moved
439  /// @param[in] iaxis the coordinate X of iatom to be moved
440  /// @return \ket{i^X} or \ket{F^\perp}
441  vecfuncT solve_cphf(const size_t iatom, const int iaxis, const Tensor<double> fock,
442  const vecfuncT& guess, const vecfuncT& rhsconst,
443  const Tensor<double> incomplete_hessian, const vecfuncT& parallel,
444  const SCFProtocol& p, const std::string& xc_data) const;
445 
446  /// solve the CPHF equation for all displacements
447 
448  /// this function computes the nemo response F^X
449  /// \f[
450  /// F^X = F^\perp + F^\parallel
451  /// \f]
452  /// To reconstruct the unregularized orbital response (not recommended):
453  /// \f[
454  /// i^X = R^X F + R F^X
455  /// \f]
456  /// The orbital response i^X constructed in this way is automatically
457  /// orthogonal to the occupied space because of the parallel term F^\parallel
458  /// @return a vector of the nemo response F^X for all displacements
459  std::vector<vecfuncT> compute_all_cphf();
460 
461  /// this function computes that part of the orbital response that is
462  /// parallel to the occupied space.
463  /// \f[
464  /// F^X = F^\perp + F^\parallel
465  /// \f]
466  /// If no NCF's are used F^\parallel vanishes.
467  /// If NCF's are used this term does not vanish because the derivatives of
468  /// the NCF does not vanish, and it is given by
469  /// \f[
470  /// F_i^\parallel = -\frac{1}{2}\sum_k|F_k ><F_k | (R^2)^X | F_i>
471  /// \f]
472  vecfuncT compute_cphf_parallel_term(const size_t iatom, const int iaxis) const;
473 
474  /// compute the IR intensities in the double harmonic approximation
475 
476  /// use the projected normal modes; units are km/mol
477  /// @param[in] normalmodes the normal modes
478  /// @param[in] dens_pt the perturbed densities for each nuclear displacement
480  const vecfuncT& dens_pt) const;
481 
482  std::shared_ptr<SCF> get_calc() const {return calc;}
483 
484  const NemoCalculationParameters& get_param() const {return param;}
485 
486  PCM get_pcm()const{return pcm;}
487 
488  /// compute the Fock matrix from scratch
489  tensorT compute_fock_matrix(const vecfuncT& nemo, const tensorT& occ) const;
490 
491  /// return a reference to the molecule
492  Molecule& molecule() {return calc->molecule;}
493 
494  /// return a reference to the molecule
495  Molecule& molecule() const {
496  return calc->molecule;
497  }
498 
499  /// make the density (alpha or beta)
501  const vecfuncT& nemo) const;
502 
503  /// make the density using different bra and ket vectors
504 
505  /// e.g. for computing the perturbed density \sum_i \phi_i \phi_i^X
506  /// or when using nemos: \sum_i R2nemo_i nemo_i
508  const vecfuncT& bra, const vecfuncT& ket, const bool refine=false) const;
509 
510  /// make the derivative of the density
511 
512  /// \f$ \nabla\rho = 2R^X R \rho_R + R^2\nabla \rho_R \f$
513  /// @param[in] rhonemo the regularized density
514  /// @param[in] axis the component of the nabla operator
515  /// @return the gradient of the *reconstructed* density
517  const int axis) const;
518 
519  /// compute the reduced densities sigma (gamma) for GGA functionals
521  const real_function_3d& rho2) const;
522 
523 
524  /// the Laplacian of the density
525 
526  /// The Laplacian should currently only be used for subsequent convolution
527  /// with a Green's function (which is reasonably stable), but not on its own!
528  ///
529  /// The Laplacian of the cuspy density is numerically fairly unstable:
530  /// - a singular term may be rewritten using the nuclear potential (see below)
531  /// - the Laplacian of the regularized density is still very noisy
532  ///
533  /// It may be computed as
534  /// \f[
535  /// \Delta \rho = \Delta (R^2 \rho_R)
536  /// = \Delta (R^2) \rho_R + 2\nabla R \nabla \rho_R + R^2 \Delta \rho_R
537  /// = 2 R^2 U1^2 \rho_R -4 R^2 ( U-V ) \rho_R + R^2 \Delta\rho_R
538  /// \f]
539  /// where we can use the identity
540  /// \f[
541  /// U=V + R^{-1}[T,R]
542  /// -2 R (U-V) = \Delta R + 2\nabla R\dot \nabla
543  /// \f]
544  /// first term comes from the definition of the U potential as the commutator
545  /// over the kinetic energy (aka the Laplacian)
546  /// @param[in] rhonemo the regularized density \rho_R
547  /// @return the laplacian of the reconstructed density \Delta (R^2\rho_R)
549 
550  /// compute the kinetic energy potential using Eq. (16) of
551  /// R. A. King and N. C. Handy, “Kinetic energy functionals from the Kohn–Sham potential,”
552  /// Phys. Chem. Chem. Phys., vol. 2, no. 22, pp. 5049–5056, 2000.
554 
555 
556  /// smooth a function by projecting it onto k-1 and then average with k
557 
558  /// kept it here for further testing
559  static void smoothen(real_function_3d& f) {
560  int k=f.get_impl()->get_k();
561  real_function_3d fproj=project(f,k-1);
562  real_function_3d freproj=project(fproj,k);
563  f=0.5*(f+freproj);
564  }
565 
566 protected:
567 
568  std::shared_ptr<SCF> calc;
569 
570 public:
572 
573 protected:
575 
576 public:
577 
578  /// return the symmetry_projector
580  return symmetry_projector;
581  }
582 
583 private:
584 
585  /// sum of square of coords at last solved geometry
586  mutable double coords_sum;
587 
588 protected:
589  /// a poisson solver
590  std::shared_ptr<real_convolution_3d> poisson;
591 
592  /// asymptotic correction for DFT
594 
595 // /// apply the AC scheme of Tozer/Handy with the multipole approximation
596 // Function<double,3> apply_ac(const Function<double,3>& vxc)const{
597 // return ac.apply(vxc);
598 // }
599 //
600 // /// apply the AC scheme of Tozer/Handy using the hartree potential
601 // Function<double,3> apply_ac(const Function<double,3>& vxc, const Function<double,3>& vhartree)const{
602 // return ac.apply(vxc,vhartree);
603 // }
604 //
605 // /// apply the AC scheme of Tozer/Handy
606 // Function<double,3> apply_ac_excited(Function<double,3>& vxc, const Function<double,3>& vhartree)const{
607 // return ac.apply_potential(vxc,vhartree);
608 // }
609 
610 private:
611  /// polarizable continuum model
613 // AC<3> ac;
614 
615 protected:
616  /// adapt the thresholds consistently to a common value
617  void set_protocol(const double thresh) {
618 
619  calc->set_protocol<3>(world,thresh);
620 
622  timer timer1(world,param.print_level()>2);
623  get_calc()->make_nuclear_potential(world);
624  construct_nuclear_correlation_factor(calc->molecule, calc->potentialmanager, param.ncf());
625  timer1.end("reproject ncf");
626  }
627 
628  // (re) construct the Poisson solver
629  poisson = std::shared_ptr<real_convolution_3d>(
631 
632  // set thresholds for the MOs
633  set_thresh(world,calc->amo,thresh);
634  set_thresh(world,calc->bmo,thresh);
635 
636  }
637 
638  /// solve the HF equations
639  double solve(const SCFProtocol& proto);
640 
641  /// given nemos, compute the HF energy using the regularized expressions for T and V
642  std::vector<double> compute_energy_regularized(const vecfuncT& nemo, const vecfuncT& Jnemo,
643  const vecfuncT& Knemo, const vecfuncT& Unemo) const;
644 
645  /// compute the reconstructed orbitals, and all potentials applied on nemo
646 
647  /// to use these potentials in the fock matrix computation they must
648  /// be multiplied by the nuclear correlation factor
649  /// @param[in] nemo the nemo orbitals
650  /// @param[out] Jnemo Coulomb operator applied on the nemos
651  /// @param[out] Knemo exchange operator applied on the nemos
652  /// @param[out] pcmnemo PCM (solvent) potential applied on the nemos
653  /// @param[out] Unemo regularized nuclear potential applied on the nemos
654  void compute_nemo_potentials(const vecfuncT& nemo,
655  vecfuncT& Jnemo, vecfuncT& Knemo, vecfuncT& xcnemo, vecfuncT& pcmnemo,
656  vecfuncT& Unemo) const;
657 
658  /// return the Coulomb potential
660 
661  /// compute the incomplete hessian
662 
663  /// incomplete hessian is the nuclear-nuclear contribution, and the
664  /// contribution from the second derivative of the nuclear potential,
665  /// and also the derivative of the nuclear correlation factor.
666  /// i.e. all contributions that *do not* contain the regularized perturbed
667  /// density, but it will contain parts of the perturbed density
669 
670  /// compute the complementary incomplete hessian
671 
672  /// @param[in] xi the response functions including the parallel part
674  const std::vector<vecfuncT>& xi) const;
675 
676  /// compute the constant term for the CPHF equations
677 
678  /// mainly all terms with the nuclear correlation factor's derivatives
679  vecfuncT make_cphf_constant_term(const size_t iatom, const int iaxis,
680  const vecfuncT& R2nemo, const real_function_3d& rhonemo) const;
681 
682 public:
683 
684  bool is_dft() const {return calc->xc.is_dft();}
685 
686  bool do_pcm() const {return param.pcm_data() != "none";}
687 
688  bool do_ac() const {return param.ac_data() != "none";}
689 
690  AC<3> get_ac() const {return ac;}
691 
692  bool do_symmetry() const {return (symmetry_projector.get_pointgroup()!="C1");}
693 
694 protected:
695 
696  /// localize the nemo orbitals
697  vecfuncT localize(const vecfuncT& nemo, const double dconv, const bool randomize) const;
698 protected:
699  /// return the threshold for vanishing elements in orbital rotations
700  double trantol() const {
701  return calc->vtol / std::min(30.0, double(get_calc()->amo.size()));
702  }
703 
704  template<typename solverT>
705  void rotate_subspace(World& world, const tensorT& U, solverT& solver,
706  int lo, int nfunc) const;
707 
708 
709  void make_plots(const real_function_3d &f,const std::string &name="function")const{
710  double width = FunctionDefaults<3>::get_cell_min_width()/2.0 - 1.e-3;
712  coord_3d start(0.0); start[0]=-width;
713  coord_3d end(0.0); end[0]=width;
714  plot_line(("line_"+name).c_str(),1000,start,end,f);
715  }
716 
717  /// save a function
718  template<typename T, size_t NDIM>
719  void save_function(const std::vector<Function<T,NDIM> >& f, const std::string name) const;
720 
721  /// load a function
722  template<typename T, size_t NDIM>
723  void load_function(std::vector<Function<T,NDIM> >& f, const std::string name) const;
724 
725 };
726 
727 /// rotate the KAIN subspace (cf. SCF.cc)
728 template<typename solverT>
729 void Nemo::rotate_subspace(World& world, const tensorT& U, solverT& solver,
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) {
734  vecfuncT& v = ulist[iter];
735  vecfuncT& r = rlist[iter];
736  vecfuncT vnew = transform(world, vecfuncT(&v[lo], &v[lo + nfunc]), U,
737  trantol(), false);
738  vecfuncT rnew = transform(world, vecfuncT(&r[lo], &r[lo + nfunc]), U,
739  trantol(), true);
740 
741  world.gop.fence();
742  for (int i=0; i<nfunc; i++) {
743  v[i] = vnew[i];
744  r[i] = rnew[i];
745  }
746  }
747  world.gop.fence();
748 }
749 
750 /// save a function
751 template<typename T, size_t NDIM>
752 void Nemo::save_function(const std::vector<Function<T,NDIM> >& f, const std::string name) const {
753  if (world.rank()==0) print("saving vector of functions",name);
755  ar & f.size();
756  for (const Function<T,NDIM>& ff:f) ar & ff;
757 }
758 
759 /// load a function
760 template<typename T, size_t NDIM>
761 void Nemo::load_function(std::vector<Function<T,NDIM> >& f, const std::string name) const {
762  if (world.rank()==0) print("loading vector of functions",name);
764  std::size_t fsize=0;
765  ar & fsize;
766  f.resize(fsize);
767  for (std::size_t i=0; i<fsize; ++i) ar & f[i];
768 }
769 
770 }
771 
772 #endif /* NEMO_H_ */
773 
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
Definition: nemo.h:69
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 &param, 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 &parallel, 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
Definition: PNO.h:27
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
Definition: TDHF.h:24
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 BinaryFstreamInputArchive.
Definition: parallel_archive.h:366
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 &param)
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.