MADNESS  0.10.1
SCF.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  $Id$
33  */
34 
35 /// \file moldft.cc
36 /// \brief Molecular HF and DFT code
37 /// \defgroup moldft The molecular density functional and Hartree-Fock code
38 
39 
40 #ifndef MADNESS_CHEM_SCF_H__INCLUDED
41 #define MADNESS_CHEM_SCF_H__INCLUDED
42 
43 #include <memory>
44 
46 #include <madness/mra/mra.h>
47 
57 #include <madness/tensor/solvers.h>
59 #include<madness/chem/pcm.h>
61 
62 #include <madness/tensor/tensor_json.hpp>
63 #include <memory>
64 
65 namespace madness {
66 
67 typedef std::shared_ptr<WorldDCPmapInterface<Key<3> > > pmapT;
69 typedef std::shared_ptr<FunctionFunctorInterface<double, 3> > functorT;
71 typedef std::vector<functionT> vecfuncT;
72 typedef std::pair<vecfuncT, vecfuncT> pairvecfuncT;
73 typedef std::vector<pairvecfuncT> subspaceT;
74 typedef Tensor<double> tensorT;
78 typedef std::shared_ptr<operatorT> poperatorT;
80 typedef std::vector<complex_functionT> cvecfuncT;
82 
83 
84 template<typename T, int NDIM>
85 struct lbcost {
86  double leaf_value;
87  double parent_value;
88 
90 
91  double operator()(const Key<NDIM>& key, const FunctionNode<T, NDIM>& node) const {
92  if (key.level() < 1) {
93  return 100.0 * (leaf_value + parent_value);
94  } else if (node.is_leaf()) {
95  return leaf_value;
96  } else {
97  return parent_value;
98  }
99  }
100 };
101 
102 
103 inline double mask1(double x) {
104  /* Iterated first beta function to switch smoothly
105  from 0->1 in [0,1]. n iterations produce 2*n-1
106  zero derivatives at the end points. Order of polyn
107  is 3^n.
108 
109  Currently use one iteration so that first deriv.
110  is zero at interior boundary and is exactly representable
111  by low order multiwavelet without refinement */
112 
113  x = (x * x * (3. - 2. * x));
114  return x;
115 }
116 
117 static double mask3(const coordT& ruser) {
118  coordT rsim;
119  user_to_sim(ruser, rsim);
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;
123 
124  if (x < lo)
125  result *= mask1(x * rlo);
126  else if (x > hi)
127  result *= mask1((1.0 - x) * rlo);
128  if (y < lo)
129  result *= mask1(y * rlo);
130  else if (y > hi)
131  result *= mask1((1.0 - y) * rlo);
132  if (z < lo)
133  result *= mask1(z * rlo);
134  else if (z > hi)
135  result *= mask1((1.0 - z) * rlo);
136 
137  return result;
138 }
139 
140 /// A MADNESS functor to compute either x, y, or z
141 class DipoleFunctor : public FunctionFunctorInterface<double, 3> {
142 private:
143  const int axis;
144 public:
146 
147  double operator()(const coordT& x) const {
148  return x[axis];
149  }
150 };
151 
152 
153 /// A MADNESS functor to compute the cartesian moment x^i * y^j * z^k (i, j, k integer and >= 0)
154 class MomentFunctor : public FunctionFunctorInterface<double, 3> {
155 private:
156  const int i, j, k;
157 public:
158  MomentFunctor(int i, int j, int k) : i(i), j(j), k(k) {}
159 
160  MomentFunctor(const std::vector<int>& x) : i(x[0]), j(x[1]), k(x[2]) {}
161 
162  double operator()(const coordT& r) const {
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];
167  return xi * yj * zk;
168  }
169 };
170 
171  class scf_data {
172 
173  std::map<std::string, std::vector<double>> e_data;
174  int iter;
175  public:
176 
177  scf_data();
178 
179  void to_json(json &j) const;
180 
181  void print_data();
182 
183  void add_data(std::map<std::string, double> values);
184  };
185 
186 
187 class SCF {
188 public:
189  std::shared_ptr<PotentialManager> potentialmanager;
190  std::shared_ptr<GTHPseudopotential<double> > gthpseudopotential;
197 
199 
200  /// alpha and beta molecular orbitals
202 
203  /// sets of orbitals grouped by their orbital energies (for localization?)
204  /// only orbitals within the same set will be mixed to localize
205  std::vector<int> aset, bset;
206 
207  /// MRA projection of the minimal basis set
209 
210  std::vector<int> at_to_bf, at_nbf;
211 
212  /// occupation numbers for alpha and beta orbitals
214 
215  /// orbital energies for alpha and beta orbitals
218  std::vector<std::shared_ptr<real_derivative_3d> > gradop;
219  double vtol;
221  double converged_for_thresh=1.e10; ///< mos are converged for this threshold
222  //double esol;//etot;
223  //double vacuo_energy;
224 
225  /// collective constructor for SCF uses contents of file \c filename and broadcasts to all nodes
226 // SCF(World & world, const char *filename);
227  /// collective constructor for SCF uses contents of stream \c input and broadcasts to all nodes
228 // SCF(World & world, std::shared_ptr<std::istream> input);
229 // SCF(World& world, const std::string& inputfile);
230  SCF(World& world, const commandlineparser& parser);
231 
232  void copy_data(World& world, const SCF& other);
233 
234  static void help() {
235  print_header2("help page for MOLDFT ");
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.");
245 
246  }
247 
248  static void print_parameters() {
250  print("default parameters for the moldft program are");
251  param.print("dft", "end");
252  print("\n\nthe molecular geometry must be specified in a separate block:");
254  }
255 
256  void set_print_timings(const bool value);
257 
258  template<std::size_t NDIM>
259  void set_protocol(World& world, double thresh) {
260  int k;
261  // Allow for imprecise conversion of threshold
262  if (thresh >= 0.9e-2)
263  k = 4;
264  else if (thresh >= 0.9e-4)
265  k = 6;
266  else if (thresh >= 0.9e-6)
267  k = 8;
268  else if (thresh >= 0.9e-8)
269  k = 10;
270  else
271  k = 12;
272 
273  // k defaults to make sense with thresh, override by providing k in input file
274  if (param.k() == -1) {
276  // param.k=k;
277  } else {
279  }
280  // don't forget to adapt the molecular smoothing parameter!! NO ... it is independent
281  // molecule.set_eprec(std::min(thresh,molecule.get_eprec()));
285  // FunctionDefaults<NDIM>::set_truncate_mode(1);
291  double safety = 0.1;
294  gradop = gradient_operator<double, 3>(world);
295 
296  // Update coefficients if using a different derivative
297  if (param.deriv() == "bspline") {
298  for (int i = 0; i < 3; ++i) (*gradop[i]).set_bspline1();
299  } else if (param.deriv() == "ble") {
300  for (int i = 0; i < 3; ++i) (*gradop[i]).set_ble1();
301  }
302 
303  mask = functionT(factoryT(world).f(mask3).initial_level(4).norefine());
304  if (world.rank() == 0 and param.print_level() > 1) {
305  print("\nSolving NDIM=", NDIM, " with thresh", thresh, " k",
307  }
308  }
309 
310  /// getter for the molecular orbitals, alpha spin
311  const vecfuncT& get_amo() const { return amo; }
312 
313  /// getter for the molecular orbitals, beta spin
314  const vecfuncT& get_bmo() const { return bmo; }
315 
316  /// getter for the occupation numbers, alpha spin
317  const tensorT& get_aocc() const { return aocc; }
318 
319  /// getter for the occupation numbers, alpha spin
320  const tensorT& get_bocc() const { return bocc; }
321 
322  bool is_spin_restricted() const { return param.get<bool>("spin_restricted"); }
323 
324  void save_mos(World& world);
325 
326  void load_mos(World& world);
327 
328  bool restart_aos(World& world);
329 
330  void do_plots(World& world);
331 
332  void project(World& world);
333 
334  void make_nuclear_potential(World& world);
335 
337 
339  const Molecule& molecule);
340 
341  void reset_aobasis(const std::string& aobasisname) {
342  aobasis = AtomicBasisSet(); // reset
343  aobasis.read_file(aobasisname);
344  }
345 
346  /// group orbitals into sets of similar orbital energies for localization
347 
348  /// @param[in] eps orbital energies
349  /// @param[in] occ occupation numbers
350  /// @param[in] nmo number of MOs for the given spin
351  /// @return vector of length nmo with the set index for each MO
352  std::vector<int> group_orbital_sets(World& world, const tensorT& eps,
353  const tensorT& occ, const int nmo) const;
354 
355  static void analyze_vectors(World& world, const vecfuncT& mo,
356  const vecfuncT& ao, double vtol,
357  const Molecule& molecule, const int print_level,
358  const AtomicBasisSet& aobasis, const tensorT& occ = tensorT(),
359  const tensorT& energy = tensorT(), const std::vector<int>& set = std::vector<int>());
360 
361  distmatT kinetic_energy_matrix(World& world, const vecfuncT& v) const;
362 
363  void initial_guess(World& world);
364 
365  void initial_load_bal(World& world);
366 
367  functionT make_density(World& world, const tensorT& occ, const vecfuncT& v) const;
368 
369  functionT make_density(World& world, const tensorT& occ, const cvecfuncT& v);
370 
371  std::vector<poperatorT> make_bsh_operators(World& world, const tensorT& evals) const;
372 
373  // Used only for initial guess that is always spin-restricted LDA
374  static functionT make_lda_potential(World& world, const functionT& arho);
375 
376 
377  // functionT make_dft_potential(World & world, const vecfuncT& vf, int ispin, int what)
378  // {
379  // return multiop_values<double, xc_potential, 3>(xc_potential(xc, ispin, what), vf);
380  // }
381 
382  double make_dft_energy(World& world, const vecfuncT& vf, int ispin) {
383  functionT vlda = multiop_values<double, xc_functional, 3>(xc_functional(xc), vf);
384  return vlda.trace();
385  }
386 
387  vecfuncT apply_potential(World& world, const tensorT& occ,
388  const vecfuncT& amo,
389  const functionT& vlocal, double& exc, double& enl, int ispin);
390 
391  tensorT derivatives(World& world, const functionT& rho) const;
392 
393  /// compute the total dipole moment of the molecule
394 
395  /// @param[in] rho the total (alpha + beta) density
396  /// @return the x,y,z components of the el. + nucl. dipole moment
397  tensorT dipole(World& world, const functionT& rho) const;
398 
399  void vector_stats(const std::vector<double>& v, double& rms,
400  double& maxabsval) const;
401 
402  vecfuncT compute_residual(World& world, tensorT& occ, tensorT& fock,
403  const vecfuncT& psi, vecfuncT& Vpsi, double& err);
404 
405  tensorT make_fock_matrix(World& world, const vecfuncT& psi,
406  const vecfuncT& Vpsi, const tensorT& occ,
407  double& ekinetic) const;
408 
409  /// make the Coulomb potential given the total density
411  return apply(*coulop, rho);
412  }
413 
414  /// Compute the two-electron integrals over the provided set of orbitals
415 
416  /// Returned is a *replicated* tensor of \f$(ij|kl)\f$ with \f$i>=j\f$
417  /// and \f$k>=l\f$. The symmetry \f$(ij|kl)=(kl|ij)\f$ is enforced.
418  Tensor<double> twoint(World& world, const vecfuncT& psi) const;
419 
420  /// compute the unitary transformation that diagonalizes the fock matrix
421 
422  /// @param[in] world the world
423  /// @param[in] overlap the overlap matrix of the orbitals
424  /// @param[in,out] fock the fock matrix; diagonal upon exit
425  /// @param[out] evals the orbital energies
426  /// @param[in] occ the occupation numbers
427  /// @param[in] thresh_degenerate threshold for orbitals being degenerate
428  /// @return the unitary matrix U: U^T F U = evals
429  tensorT get_fock_transformation(World& world, const tensorT& overlap,
430  tensorT& fock, tensorT& evals, const tensorT& occ,
431  const double thresh_degenerate) const;
432 
433 
434  /// diagonalize the fock matrix, taking care of degenerate states
435 
436  /// Vpsi is passed in to make sure orbitals and Vpsi are in phase
437  /// @param[in] world the world
438  /// @param[in,out] fock the fock matrix (diagonal upon exit)
439  /// @param[in,out] psi the orbitals
440  /// @param[in,out] Vpsi the orbital times the potential
441  /// @param[out] evals the orbital energies
442  /// @param[in] occ occupation numbers
443  /// @param[in] thresh threshold for rotation and truncation
444  /// @return the unitary matrix U: U^T F U = evals
445  tensorT diag_fock_matrix(World& world, tensorT& fock,
446  vecfuncT& psi, vecfuncT& Vpsi, tensorT& evals,
447  const tensorT& occ, const double thresh) const;
448 
449 
450  void loadbal(World& world, functionT& arho, functionT& brho, functionT& arho_old,
451  functionT& brho_old, subspaceT& subspace);
452 
453 
454  void rotate_subspace(World& world, const tensorT& U, subspaceT& subspace,
455  int lo, int nfunc, double trantol) const;
456 
457  void rotate_subspace(World& world, const distmatT& U, subspaceT& subspace,
458  int lo, int nfunc, double trantol) const;
459 
460  void update_subspace(World& world,
461  vecfuncT& Vpsia, vecfuncT& Vpsib,
462  tensorT& focka, tensorT& fockb,
464  double& bsh_residual, double& update_residual);
465 
466  /// perform step restriction following the KAIN solver
467 
468  /// undo the rotation from the KAIN solver if the rotation exceeds the
469  /// maxrotn parameter
470  /// @param[in] world the world
471  /// @param[in] mo vector of orbitals from previous iteration
472  /// @param[in,out] mo_new vector of orbitals from the KAIN solver
473  /// @param[in] spin "alpha" or "beta" for user information
474  /// @return max residual
475  double do_step_restriction(World& world, const vecfuncT& mo,
476  vecfuncT& mo_new, std::string spin) const;
477 
478  /// orthonormalize the vectors
479 
480  /// @param[in] world the world
481  /// @param[in,out] amo_new the vectors to be orthonormalized
482  void orthonormalize(World& world, vecfuncT& amo_new) const;
483 
484  void orthonormalize(World& world, vecfuncT& amo_new, int nocc) const;
485 
487  complex_functionT r = psi; // Shallow copy violates constness !!!!!!!!!!!!!!!!!
488  coordT lo, hi;
489  lo[2] = -10;
490  hi[2] = +10;
491 
492  r.reconstruct();
493  r.broaden();
494  r.broaden();
495  r.broaden();
496  r.broaden();
497  r = apply_1d_realspace_push(*q1d, r, 2);
498  r.sum_down();
499  r = apply_1d_realspace_push(*q1d, r, 1);
500  r.sum_down();
501  r = apply_1d_realspace_push(*q1d, r, 0);
502  r.sum_down();
503 
504  return r;
505  }
506 
507  // For given protocol, solve the DFT/HF/response equations
508  void solve(World& world);
509 
510  void output_calc_info_schema() const;
511 
512 // void output_scf_info_schema(const std::map<std::string, double> &vals,
513 // const tensorT &dipole_T) const;
514 
515 };
516 
517 // Computes molecular energy as a function of the geometry
518 // This is cludgy ... need better factorization of functionality
519 // between calculation, main program and this ... or just merge it all.
523  mutable double coords_sum; // sum of square of coords at last solved geometry
524 
525 public:
527  : world(world), calc(calc), coords_sum(-1.0) {}
528 
529  std::string name() const { return "Molecularenerg"; }
530 
531  bool selftest() { return true; }
532 
533  bool provides_gradient() const { return true; }
534 
535  double value(const Tensor<double>& x) {
536  double xsq = x.sumsq();
537  if (xsq == coords_sum) {
538  return calc.current_energy;
539  }
541  coords_sum = xsq;
542 
543  // read converged wave function from disk if there is one
544  if (calc.param.no_compute()) {
548  return calc.current_energy;
549  }
550 
551  // initialize the PCM solver for this geometry
552  if (calc.param.pcm_data() != "none") {
554  }
555 
556  bool have_initial_guess = false;
557  // get initial orbitals
558  if (calc.param.restart()) {
560  have_initial_guess = true;
561  } else if (calc.param.restartao()) {
562  calc.reset_aobasis("sto-3g");
564  have_initial_guess = calc.restart_aos(world);
565  }
566 
567  if (not have_initial_guess) {
572  }
573 
574  calc.reset_aobasis("sto-3g");
575  calc.ao.clear(); world.gop.fence();
577 
578 
579  // The below is missing convergence test logic, etc.
580 
581  // Make the nuclear potential, initial orbitals, etc.
582  for (unsigned int proto = 0; proto < calc.param.protocol().size(); proto++) {
583 
584  int nvalpha = calc.param.nmo_alpha() - calc.param.nalpha();
585  int nvbeta = calc.param.nmo_beta() - calc.param.nbeta();
586  int nvalpha_start, nv_old;
587 
588  //repeat with gradually decreasing nvirt, only for first protocol
589  if (proto == 0 && nvalpha > 0) {
590  nvalpha_start = nvalpha * calc.param.nv_factor();
591  } else {
592  nvalpha_start = nvalpha;
593  }
594 
595  nv_old = nvalpha_start;
596 
597  for (int nv = nvalpha_start; nv >= nvalpha; nv -= nvalpha) {
598 
599  if (nv > 0 && world.rank() == 0) std::cout << "Running with " << nv << " virtual states" << std::endl;
600 
601  calc.param.set_user_defined_value("nmo_alpha", calc.param.nalpha() + nv);
602  // check whether this is sensible for spin restricted case
603  if (calc.param.nbeta() && !calc.param.spin_restricted()) {
604  if (nvbeta == nvalpha) {
605  calc.param.set_user_defined_value("nmo_beta", calc.param.nbeta() + nv);
606  } else {
607  calc.param.set_user_defined_value("nmo_beta", calc.param.nbeta() + nv + nvbeta - nvalpha);
608  }
609  }
610 
611  calc.set_protocol<3>(world, calc.param.protocol()[proto]);
613 
614  if (nv != nv_old) {
615  calc.amo.resize(calc.param.nmo_alpha());
616  calc.bmo.resize(calc.param.nmo_beta());
617 
619  for (int i = 0; i < calc.param.nalpha(); ++i)
620  calc.aocc[i] = 1.0;
621 
623  for (int i = 0; i < calc.param.nbeta(); ++i)
624  calc.bocc[i] = 1.0;
625 
626  // might need to resize aset, bset, but for the moment this doesn't seem to be necessary
627 
628  }
629 
630  // project orbitals into higher k
631  if (proto > 0) calc.project(world);
632 
633  // If the basis for the inital guess was not sto-3g
634  // switch to sto-3g since this is needed for analysis
635  // of the MOs and orbital localization
636  // Only do this if not starting from NWChem.
637  // analysis will be done on NWChem orbitals.
638 
639  if (calc.param.aobasis() != "sto-3g") { // was also && calc.param.nwfile() == "none"
640  calc.reset_aobasis("sto-3g");
641  }
642  calc.ao.clear(); world.gop.fence();
644  calc.solve(world);
645 
646  if (calc.param.save())
648 
649  nv_old = nv;
650  // exit loop over decreasing nvirt if nvirt=0
651  if (nv == 0) break;
652 
653  }
654 
655  }
656  return calc.current_energy;
657  }
658 
660  value(x); // Ensures DFT equations are solved at this geometry
661 
663  functionT brho = rho;
664  if (!calc.param.spin_restricted())
665  brho = calc.make_density(world, calc.bocc, calc.bmo);
666  rho.gaxpy(1.0, brho, 1.0);
667 
668  return calc.derivatives(world, rho);
669  }
670 
671 
673  value(molecule.get_all_coords().flat()); // Ensures DFT equations are solved at this geometry
674 
676  functionT brho = rho;
677  if (!calc.param.spin_restricted())
678  brho = calc.make_density(world, calc.bocc, calc.bmo);
679  rho.gaxpy(1.0, brho, 1.0);
680 
682  gradient = calc.derivatives(world, rho);
683  }
684 
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;
690 
692 
693  nlohmann::json calc_precision={ };
694  calc_precision["eprec"]=calc.molecule.parameters.eprec();
695  calc_precision["dconv"]=calc.param.dconv();
696  calc_precision["econv"]=calc.param.econv();
697  calc_precision["thresh"]=FunctionDefaults<3>::get_thresh();
698  calc_precision["k"]=FunctionDefaults<3>::get_k();
699 
700  auto mol_json=this->calc.molecule.to_json();
701 
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()});
705  int_vals.push_back({"calcinfo_natom", calc.molecule.natom()});
706 
707 
708  to_json(j, int_vals);
709  double_vals.push_back({"return_energy", value(calc.molecule.get_all_coords().flat())});
710  to_json(j, double_vals);
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});
714  }
715 
716  to_json(j, double_tensor_vals);
717  param.to_json(j);
718  calc.e_data.to_json(j);
719 
720  j["precision"]=calc_precision;
721  j["molecule"]=mol_json;
722 
723  output_schema(param.prefix()+".calc_info", j);
724  }
725 
726 
727 };
728 }
729 
730 #endif /* SCF_H_ */
731 
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
FunctionFunctorInterface()
Definition: function_interface.h:77
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:209
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
Definition: SCF.h:520
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
madness::Tensor< double > gradient(const Tensor< double > &x)
Should return the derivative of the function.
Definition: SCF.h:659
double coords_sum
Definition: SCF.h:523
void output_calc_info_schema()
Definition: SCF.h:685
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
Definition: SCF.h:187
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
const tensorT & get_aocc() const
getter for the occupation numbers, alpha spin
Definition: SCF.h:317
static vecfuncT project_ao_basis_only(World &world, const AtomicBasisSet &aobasis, const Molecule &molecule)
Definition: SCF.cc:581
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
SCF(World &world, const commandlineparser &parser)
collective constructor for SCF uses contents of file filename and broadcasts to all nodes
Definition: SCF.cc:234
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 vecfuncT & get_amo() const
getter for the molecular orbitals, alpha spin
Definition: SCF.h:311
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
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
const vecfuncT & get_bmo() const
getter for the molecular orbitals, beta spin
Definition: SCF.h:314
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
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
const tensorT & get_bocc() const
getter for the occupation numbers, alpha spin
Definition: SCF.h:320
Tensor< T > reshape(int ndimnew, const long *d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition: tensor.h:1384
Tensor< T > flat()
Returns new view/tensor rehshaping to flat (1-d) tensor.
Definition: tensor.h:1555
T sumsq() const
Returns the sum of the squares of the elements.
Definition: tensor.h:1669
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
Definition: SCF.h:171
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
#define max(a, b)
Definition: lda.h:51
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
void print_header2(const std::string &s)
medium section heading
Definition: print.cc:54
Vector< double, 3 > coordT
Definition: corepotential.cc:54
nlohmann::json json
Definition: QCCalculationParametersBase.h:27
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
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
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
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
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
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
double Q(double a)
Definition: relops.cc:20
nlohmann::json json
Definition: response_parameters.cpp:6
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.
std::string prefix
Definition: tdse.cc:71
Definition: CalculationParameters.h:51
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::vector< double > protocol() const
Definition: CalculationParameters.h:202
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
Definition: SCF.h:85
double parent_value
Definition: SCF.h:87
double leaf_value
Definition: SCF.h:86
lbcost(double leaf_value=1.0, double parent_value=0.0)
Definition: SCF.h:89
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