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
59#include<madness/chem/pcm.h>
61
62#include <madness/tensor/tensor_json.hpp>
63#include <memory>
64
65namespace madness {
66
67typedef std::shared_ptr<WorldDCPmapInterface<Key<3> > > pmapT;
69typedef std::shared_ptr<FunctionFunctorInterface<double, 3> > functorT;
71typedef std::vector<functionT> vecfuncT;
72typedef std::pair<vecfuncT, vecfuncT> pairvecfuncT;
73typedef std::vector<pairvecfuncT> subspaceT;
78typedef std::shared_ptr<operatorT> poperatorT;
80typedef std::vector<complex_functionT> cvecfuncT;
82
83
84template<typename T, int NDIM>
85struct lbcost {
86 double leaf_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
103inline 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
117static double mask3(const coordT& ruser) {
118 coordT 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
141class DipoleFunctor : public FunctionFunctorInterface<double, 3> {
142private:
143 const int axis;
144public:
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)
154class MomentFunctor : public FunctionFunctorInterface<double, 3> {
155private:
156 const int i, j, k;
157public:
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
187class SCF {
188public:
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;
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",
306 FunctionDefaults<NDIM>::get_k(), " conv", std::max(thresh, param.dconv()), "\n");
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
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 get_initial_orbitals(World& world);
364
365 void initial_guess(World& world);
366
367 void initial_guess_from_nwchem(World& world);
368
369 void initial_load_bal(World& world);
370
371 functionT make_density(World& world, const tensorT& occ, const vecfuncT& v) const;
372
373 functionT make_density(World& world, const tensorT& occ, const cvecfuncT& v);
374
375 static std::vector<poperatorT> make_bsh_operators(World& world, const tensorT& evals,
377
378 // Used only for initial guess that is always spin-restricted LDA
379 static functionT make_lda_potential(World& world, const functionT& arho);
380
381
382 // functionT make_dft_potential(World & world, const vecfuncT& vf, int ispin, int what)
383 // {
384 // return multiop_values<double, xc_potential, 3>(xc_potential(xc, ispin, what), vf);
385 // }
386
387 double make_dft_energy(World& world, const vecfuncT& vf, int ispin) {
389 return vlda.trace();
390 }
391
392 vecfuncT apply_potential(World& world, const tensorT& occ,
393 const vecfuncT& amo,
394 const functionT& vlocal, double& exc, double& enl, int ispin);
395
396 tensorT derivatives(World& world, const functionT& rho) const;
397
398 /// compute the total dipole moment of the molecule
399
400 /// @param[in] rho the total (alpha + beta) density
401 /// @return the x,y,z components of the el. + nucl. dipole moment
402 tensorT dipole(World& world, const functionT& rho) const;
403
404 void vector_stats(const std::vector<double>& v, double& rms,
405 double& maxabsval) const;
406
407 vecfuncT compute_residual(World& world, tensorT& occ, tensorT& fock,
408 const vecfuncT& psi, vecfuncT& Vpsi, double& err);
409
411 const vecfuncT& Vpsi, const tensorT& occ,
412 double& ekinetic) const;
413
414 /// make the Coulomb potential given the total density
416 return apply(*coulop, rho);
417 }
418
419 /// Compute the two-electron integrals over the provided set of orbitals
420
421 /// Returned is a *replicated* tensor of \f$(ij|kl)\f$ with \f$i>=j\f$
422 /// and \f$k>=l\f$. The symmetry \f$(ij|kl)=(kl|ij)\f$ is enforced.
423 Tensor<double> twoint(World& world, const vecfuncT& psi) const;
424
425 /// compute the unitary transformation that diagonalizes the fock matrix
426
427 /// @param[in] world the world
428 /// @param[in] overlap the overlap matrix of the orbitals
429 /// @param[in,out] fock the fock matrix; diagonal upon exit
430 /// @param[out] evals the orbital energies
431 /// @param[in] occ the occupation numbers
432 /// @param[in] thresh_degenerate threshold for orbitals being degenerate
433 /// @return the unitary matrix U: U^T F U = evals
435 tensorT& fock, tensorT& evals, const tensorT& occ,
436 const double thresh_degenerate) const;
437
438
439 /// diagonalize the fock matrix, taking care of degenerate states
440
441 /// Vpsi is passed in to make sure orbitals and Vpsi are in phase
442 /// @param[in] world the world
443 /// @param[in,out] fock the fock matrix (diagonal upon exit)
444 /// @param[in,out] psi the orbitals
445 /// @param[in,out] Vpsi the orbital times the potential
446 /// @param[out] evals the orbital energies
447 /// @param[in] occ occupation numbers
448 /// @param[in] thresh threshold for rotation and truncation
449 /// @return the unitary matrix U: U^T F U = evals
452 const tensorT& occ, const double thresh) const;
453
454
457
458
459 void rotate_subspace(World& world, const tensorT& U, subspaceT& subspace,
460 int lo, int nfunc, double trantol) const;
461
462 void rotate_subspace(World& world, const distmatT& U, subspaceT& subspace,
463 int lo, int nfunc, double trantol) const;
464
465 void update_subspace(World& world,
469 double& bsh_residual, double& update_residual);
470
471 /// perform step restriction following the KAIN solver
472
473 /// undo the rotation from the KAIN solver if the rotation exceeds the
474 /// maxrotn parameter
475 /// @param[in] world the world
476 /// @param[in] mo vector of orbitals from previous iteration
477 /// @param[in,out] mo_new vector of orbitals from the KAIN solver
478 /// @param[in] spin "alpha" or "beta" for user information
479 /// @return max residual
480 double do_step_restriction(World& world, const vecfuncT& mo,
481 vecfuncT& mo_new, std::string spin) const;
482
483 /// orthonormalize the vectors
484
485 /// @param[in] world the world
486 /// @param[in,out] amo_new the vectors to be orthonormalized
487 void orthonormalize(World& world, vecfuncT& amo_new) const;
488
489 void orthonormalize(World& world, vecfuncT& amo_new, int nocc) const;
490
492 complex_functionT r = psi; // Shallow copy violates constness !!!!!!!!!!!!!!!!!
493 coordT lo, hi;
494 lo[2] = -10;
495 hi[2] = +10;
496
497 r.reconstruct();
498 r.broaden();
499 r.broaden();
500 r.broaden();
501 r.broaden();
502 r = apply_1d_realspace_push(*q1d, r, 2);
503 r.sum_down();
504 r = apply_1d_realspace_push(*q1d, r, 1);
505 r.sum_down();
506 r = apply_1d_realspace_push(*q1d, r, 0);
507 r.sum_down();
508
509 return r;
510 }
511
512 // For given protocol, solve the DFT/HF/response equations
513 void solve(World& world);
514
515 void output_calc_info_schema() const;
516
517// void output_scf_info_schema(const std::map<std::string, double> &vals,
518// const tensorT &dipole_T) const;
519
520};
521
522// Computes molecular energy as a function of the geometry
523// This is cludgy ... need better factorization of functionality
524// between calculation, main program and this ... or just merge it all.
528 mutable double coords_sum; // sum of square of coords at last solved geometry
529
530public:
533
534 std::string name() const { return "Molecularenerg"; }
535
536 bool selftest() { return true; }
537
538 bool provides_gradient() const { return true; }
539
540 double value(const Tensor<double>& x) {
541 double xsq = x.sumsq();
542 if (xsq == coords_sum) {
543 return calc.current_energy;
544 }
546 coords_sum = xsq;
547
548 // read converged wave function from disk if there is one
549 if (calc.param.no_compute()) {
553 return calc.current_energy;
554 }
555
556 // initialize the PCM solver for this geometry
557 if (calc.param.pcm_data() != "none") {
559 }
560
562
563 // AOs are needed for final analysis, and for localization
564 calc.reset_aobasis("sto-3g");
565 calc.ao.clear(); world.gop.fence();
567
568 // The below is missing convergence test logic, etc.
569
570 // Make the nuclear potential, initial orbitals, etc.
571 for (unsigned int proto = 0; proto < calc.param.protocol().size(); proto++) {
572
573 int nvalpha = calc.param.nmo_alpha() - calc.param.nalpha();
574 int nvbeta = calc.param.nmo_beta() - calc.param.nbeta();
576
577 //repeat with gradually decreasing nvirt, only for first protocol
578 if (proto == 0 && nvalpha > 0) {
579 nvalpha_start = nvalpha * calc.param.nv_factor();
580 } else {
581 nvalpha_start = nvalpha;
582 }
583
585
586 for (int nv = nvalpha_start; nv >= nvalpha; nv -= nvalpha) {
587
588 if (nv > 0 && world.rank() == 0) std::cout << "Running with " << nv << " virtual states" << std::endl;
589
591 // check whether this is sensible for spin restricted case
593 if (nvbeta == nvalpha) {
595 } else {
596 calc.param.set_user_defined_value("nmo_beta", calc.param.nbeta() + nv + nvbeta - nvalpha);
597 }
598 }
599
602
603 if (nv != nv_old) {
604 calc.amo.resize(calc.param.nmo_alpha());
605 calc.bmo.resize(calc.param.nmo_beta());
606
608 for (int i = 0; i < calc.param.nalpha(); ++i)
609 calc.aocc[i] = 1.0;
610
612 for (int i = 0; i < calc.param.nbeta(); ++i)
613 calc.bocc[i] = 1.0;
614
615 // might need to resize aset, bset, but for the moment this doesn't seem to be necessary
616
617 }
618
619 // project orbitals into higher k
620 if (proto > 0) calc.project(world);
621
622 // If the basis for the inital guess was not sto-3g
623 // switch to sto-3g since this is needed for analysis
624 // of the MOs and orbital localization
625 // Only do this if not starting from NWChem.
626 // analysis will be done on NWChem orbitals.
627
628 if (calc.param.aobasis() != "sto-3g") { // was also && calc.param.nwfile() == "none"
629 calc.reset_aobasis("sto-3g");
630 }
631 calc.ao.clear(); world.gop.fence();
634
635 if (calc.param.save())
637
638 nv_old = nv;
639 // exit loop over decreasing nvirt if nvirt=0
640 if (nv == 0) break;
641
642 }
643
644 }
645 return calc.current_energy;
646 }
647
649 value(x); // Ensures DFT equations are solved at this geometry
650
652 functionT brho = rho;
655 rho.gaxpy(1.0, brho, 1.0);
656
657 return calc.derivatives(world, rho);
658 }
659
660
662 value(molecule.get_all_coords().flat()); // Ensures DFT equations are solved at this geometry
663
665 functionT brho = rho;
668 rho.gaxpy(1.0, brho, 1.0);
669
672 }
673
675 nlohmann::json j = {};
679
681
682 nlohmann::json calc_precision={ };
684 calc_precision["dconv"]=calc.param.dconv();
685 calc_precision["econv"]=calc.param.econv();
688
689 auto mol_json=this->calc.molecule.to_json();
690
691 int_vals.push_back({"calcinfo_nmo", param.nmo_alpha() + param.nmo_beta()});
692 int_vals.push_back({"calcinfo_nalpha", param.nalpha()});
693 int_vals.push_back({"calcinfo_nbeta", param.nbeta()});
694 int_vals.push_back({"calcinfo_natom", calc.molecule.natom()});
695
696
697 to_json(j, int_vals);
698 double_vals.push_back({"return_energy", value(calc.molecule.get_all_coords().flat())});
700 double_tensor_vals.push_back({"scf_eigenvalues_a", calc.aeps});
701 if (param.nbeta() != 0 && !param.spin_restricted()) {
702 double_tensor_vals.push_back({"scf_eigenvalues_b", calc.beps});
703 }
704
706 param.to_json(j);
708
709 j["precision"]=calc_precision;
710 j["molecule"]=mol_json;
711
712 output_schema(param.prefix()+".calc_info", j);
713 }
714
715
716};
717}
718
719#endif /* SCF_H_ */
720
Operators for the molecular HF and DFT code.
A MADNESS functor to compute either x, y, or z.
Definition preal.cc:115
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:258
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:100
static int get_k()
Returns the default wavelet order.
Definition funcdefaults.h:163
static void set_apply_randomize(bool value)
Sets the random load balancing for integral operators flag.
Definition funcdefaults.h:294
static void set_thresh(double value)
Sets the default threshold.
Definition funcdefaults.h:183
static void set_k(int value)
Sets the default wavelet order.
Definition funcdefaults.h:170
static const double & get_thresh()
Returns the default threshold.
Definition funcdefaults.h:176
static void set_autorefine(bool value)
Sets the default adaptive autorefinement flag.
Definition funcdefaults.h:260
static void set_project_randomize(bool value)
Sets the random load balancing for projection flag.
Definition funcdefaults.h:305
static void set_initial_level(int value)
Sets the default initial projection level.
Definition funcdefaults.h:200
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:362
static void set_refine(bool value)
Sets the default adaptive refinement flag.
Definition funcdefaults.h:248
FunctionFactory implements the named-parameter idiom for Function.
Definition function_factory.h:86
Abstract base class interface required for functors used as input to Functions.
Definition function_interface.h:68
FunctionNode holds the coefficients, etc., at each node of the 2^NDIM-tree.
Definition funcimpl.h:127
A multiresolution adaptive numerical function.
Definition mra.h:139
void broaden(const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), bool fence=true) const
Inplace broadens support in scaling function basis.
Definition mra.h:861
void sum_down(bool fence=true) const
Sums scaling coeffs down tree restoring state with coeffs only at leaves. Optional fence....
Definition mra.h:823
const Function< T, NDIM > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition mra.h:800
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:1006
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:68
Level level() const
Definition key.h:161
Definition SCF.h:525
bool selftest()
Definition SCF.h:536
bool provides_gradient() const
Override this to return true if the derivative is implemented.
Definition SCF.h:538
double value(const Tensor< double > &x)
Should return the value of the objective function.
Definition SCF.h:540
World & world
Definition SCF.h:526
SCF & calc
Definition SCF.h:527
double coords_sum
Definition SCF.h:528
void output_calc_info_schema()
Definition SCF.h:674
madness::Tensor< double > gradient(const Tensor< double > &x)
Should return the derivative of the function.
Definition SCF.h:648
void energy_and_gradient(const Molecule &molecule, double &energy, Tensor< double > &gradient)
Definition SCF.h:661
std::string name() const
Definition SCF.h:534
MolecularEnergy(World &world, SCF &calc)
Definition SCF.h:531
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:545
tensorT derivatives(World &world, const functionT &rho) const
Definition SCF.cc:1464
void save_mos(World &world)
Definition SCF.cc:318
static vecfuncT project_ao_basis_only(World &world, const AtomicBasisSet &aobasis, const Molecule &molecule)
Definition SCF.cc:639
const tensorT & get_aocc() const
getter for the occupation numbers, alpha spin
Definition SCF.h:317
std::shared_ptr< GTHPseudopotential< double > > gthpseudopotential
Definition SCF.h:190
void make_nuclear_potential(World &world)
Definition SCF.cc:614
void get_initial_orbitals(World &world)
get the initial orbitals for a calculation
Definition SCF.cc:495
std::vector< int > aset
Definition SCF.h:205
static void print_parameters()
Definition SCF.h:248
AtomicBasisSet aobasis
Definition SCF.h:195
vecfuncT ao
MRA projection of the minimal basis set.
Definition SCF.h:208
scf_data e_data
Definition SCF.h:198
void initial_guess(World &world)
Definition SCF.cc:1027
vecfuncT apply_potential(World &world, const tensorT &occ, const vecfuncT &amo, const functionT &vlocal, double &exc, double &enl, int ispin)
Definition SCF.cc:1384
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:1575
const tensorT & get_bocc() const
getter for the occupation numbers, alpha spin
Definition SCF.h:320
std::vector< int > at_to_bf
Definition SCF.h:210
std::vector< int > bset
Definition SCF.h:205
distmatT kinetic_energy_matrix(World &world, const vecfuncT &v) const
Definition SCF.cc:711
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:1537
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:387
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:1892
Molecule molecule
Definition SCF.h:191
void load_mos(World &world)
Definition SCF.cc:358
tensorT beps
Definition SCF.h:216
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:1790
const vecfuncT & get_amo() const
getter for the molecular orbitals, alpha spin
Definition SCF.h:311
static void analyze_vectors(World &world, const vecfuncT &mo, const vecfuncT &ao, double vtol, const Molecule &molecule, const int print_level, const AtomicBasisSet &aobasis, const tensorT &occ=tensorT(), const tensorT &energy=tensorT(), const std::vector< int > &set=std::vector< int >())
Definition SCF.cc:653
XCfunctional xc
Definition SCF.h:193
void initial_load_bal(World &world)
Definition SCF.cc:1294
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:2034
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:1759
bool restart_aos(World &world)
Definition SCF.cc:741
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:415
tensorT make_fock_matrix(World &world, const vecfuncT &psi, const vecfuncT &Vpsi, const tensorT &occ, double &ekinetic) const
Definition SCF.cc:1671
Tensor< double > twoint(World &world, const vecfuncT &psi) const
Compute the two-electron integrals over the provided set of orbitals.
Definition SCF.cc:1729
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:1376
functionT make_density(World &world, const tensorT &occ, const vecfuncT &v) const
Definition SCF.cc:1311
void rotate_subspace(World &world, const tensorT &U, subspaceT &subspace, int lo, int nfunc, double trantol) const
Definition SCF.cc:1858
void reset_aobasis(const std::string &aobasisname)
Definition SCF.h:341
void project(World &world)
Definition SCF.cc:592
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:1823
void vector_stats(const std::vector< double > &v, double &rms, double &maxabsval) const
Definition SCF.cc:1563
vecfuncT project_ao_basis(World &world, const AtomicBasisSet &aobasis)
Definition SCF.cc:631
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:1269
complex_functionT APPLY(const complex_operatorT *q1d, const complex_functionT &psi)
Definition SCF.h:491
void initial_guess_from_nwchem(World &world)
Definition SCF.cc:794
functionT mask
Definition SCF.h:196
std::shared_ptr< PotentialManager > potentialmanager
Definition SCF.h:189
static std::vector< poperatorT > make_bsh_operators(World &world, const tensorT &evals, const CalculationParameters &param)
Definition SCF.cc:1352
const vecfuncT & get_bmo() const
getter for the molecular orbitals, beta spin
Definition SCF.h:314
vecfuncT bmo
Definition SCF.h:201
void solve(World &world)
Definition SCF.cc:2156
void orthonormalize(World &world, vecfuncT &amo_new) const
orthonormalize the vectors
Definition SCF.cc:2106
Convolutions in separated form (including Gaussian)
Definition operator.h:139
A tensor is a multidimensional array.
Definition tensor.h:317
Tensor< T > reshape(int ndimnew, const long *d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition tensor.h:1384
T sumsq() const
Returns the sum of the squares of the elements.
Definition tensor.h:1669
Tensor< T > flat()
Returns new view/tensor rehshaping to flat (1-d) tensor.
Definition tensor.h:1555
A simple, fixed dimension vector.
Definition vector.h:64
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition worldgop.cc:161
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:320
WorldGopInterface & gop
Global operations.
Definition world.h:207
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
Main include file for MADNESS and defines Function interface.
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
void print_header2(const std::string &s)
medium section heading
Definition print.cc:54
Function< TENSOR_RESULT_TYPE(typename opT::opT, R), NDIM > apply_1d_realspace_push(const opT &op, const Function< R, NDIM > &f, int axis, bool fence=true)
Definition mra.h:2263
Vector< double, 3 > coordT
Definition corepotential.cc:54
static SeparatedConvolution< double, 3 > * CoulombOperatorPtr(World &world, double lo, double eps, const array_of_bools< 3 > &lattice_sum=FunctionDefaults< 3 >::get_bc().is_periodic(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition operator.h:1818
nlohmann::json json
Definition QCCalculationParametersBase.h:27
static void user_to_sim(const Vector< double, NDIM > &xuser, Vector< double, NDIM > &xsim)
Convert user coords (cell[][]) to simulation coords ([0,1]^ndim)
Definition funcdefaults.h:425
Tensor< double > tensorT
Definition distpm.cc:21
Function< std::complex< double >, 3 > complex_functionT
Definition SCF.h:79
std::vector< pairvecfuncT > subspaceT
Definition SCF.h:73
double mask1(double x)
Definition SCF.h:103
DistributedMatrix< double > distmatT
Definition SCF.h:75
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition SCF.h:72
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition print.h:225
FunctionFactory< double, 3 > factoryT
Definition corepotential.cc:57
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
std::shared_ptr< operatorT> poperatorT
Definition SCF.h:78
Function< double, 3 > functionT
Definition corepotential.cc:56
NDIM & f
Definition mra.h:2448
void to_json(nlohmann::json &j)
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition corepotential.cc:55
std::vector< complex_functionT > cvecfuncT
Definition SCF.h:80
vector< functionT > vecfuncT
Definition corepotential.cc:58
static double mask3(const coordT &ruser)
Definition SCF.h:117
std::shared_ptr< WorldDCPmapInterface< Key< 3 > > > pmapT
Definition SCF.h:67
Convolution1D< double_complex > complex_operatorT
Definition SCF.h:81
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
SeparatedConvolution< double, 3 > operatorT
Definition SCF.h:77
double Q(double a)
Definition relops.cc:20
static const double thresh
Definition rk.cc:45
static const long k
Definition rk.cc:44
const double xi
Exponent for delta function approx.
Definition siam_example.cc:60
Defines interfaces for optimization and non-linear equation solvers.
std::string prefix
Definition tdse.cc:71
Definition CalculationParameters.h:51
std::vector< double > protocol() const
Definition CalculationParameters.h:202
int nv_factor() const
Definition CalculationParameters.h:167
int k() const
Definition CalculationParameters.h:179
int nmo_alpha() const
Definition CalculationParameters.h:169
bool save() const
Definition CalculationParameters.h:203
double dconv() const
Definition CalculationParameters.h:142
double econv() const
Definition CalculationParameters.h:141
int print_level() const
Definition CalculationParameters.h:188
std::string deriv() const
Definition CalculationParameters.h:193
double L() const
Definition CalculationParameters.h:178
std::string aobasis() const
Definition CalculationParameters.h:200
int nalpha() const
Definition CalculationParameters.h:162
int nbeta() const
Definition CalculationParameters.h:163
int nmo_beta() const
Definition CalculationParameters.h:170
bool no_compute() const
Definition CalculationParameters.h:175
bool spin_restricted() const
Definition CalculationParameters.h:174
double lo() const
Definition CalculationParameters.h:177
std::string pcm_data() const
Definition CalculationParameters.h:195
Definition convolution1d.h:990
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
lbcost(double leaf_value=1.0, double parent_value=0.0)
Definition SCF.h:89
double leaf_value
Definition SCF.h:86
double operator()(const Key< NDIM > &key, const FunctionNode< T, NDIM > &node) const
Definition SCF.h:91
Class to compute the energy functional.
Definition xcfunctional.h:360
InputParameters param
Definition tdse.cc:203
constexpr std::size_t NDIM
Definition testgconv.cc:54
static Molecule molecule
Definition testperiodicdft.cc:39
static Subspace * subspace
Definition testperiodicdft.cc:41