MADNESS 0.10.1
MolecularOrbitals.h
Go to the documentation of this file.
1/*
2 * MolecularOrbitals.h
3 *
4 * Created on: 11 Jul 2019
5 * Author: fbischoff
6 */
7
8#ifndef SRC_APPS_CHEM_MOLECULARORBITALS_H_
9#define SRC_APPS_CHEM_MOLECULARORBITALS_H_
10
11#include<vector>
15#include<stdio.h>
16
17namespace madness {
18
19template<typename T, std::size_t NDIM> class Function;
20//class World;
21class AtomicBasisSet;
22
23template<typename T, std::size_t NDIM>
25public:
26
28
29 MolecularOrbitals() = default;
30
32 : mo(mo), eps(), irreps(), occ(), localize_sets() {
33 }
34
36 : mo(mo), eps(eps), irreps(), occ(), localize_sets() {
37 }
38
40 const std::vector<std::string>& irrep, const Tensor<double>& occ,
41 const std::vector<int>& set)
42 : mo(mo), eps(eps), irreps(irrep), occ(occ), localize_sets(set) {
43 }
44
45 MolecularOrbitals get_subset(const int iset) const {
47 auto s=slices[iset];
48 MolecularOrbitals result;
49 MADNESS_CHECK(mo.size()>=size_t(s.end+1));
50// result.mo.assign(mo.begin()+s.start,mo.begin()+s.end+1);
51 for (int i=s.start; i<s.end+1; ++i) result.mo.push_back(copy(mo[i]));
52 if (eps.size()>0) result.eps=copy(eps(s));
53 if (irreps.size()>0) result.irreps.assign(irreps.begin()+s.start,irreps.begin()+s.end+1);
54 if (occ.size()>0) result.occ=copy(occ(s));
55 if (localize_sets.size()>0) result.localize_sets.assign(localize_sets.begin()+s.start,localize_sets.begin()+s.end+1);
56
57 return result;
58 }
59
60 std::size_t size() const {
61 return mo.size();
62 }
63
64 std::vector<Function<T,NDIM> > get_mos() const {
65 return mo;
66 }
67
68 [[nodiscard]] Tensor<double> get_eps() const {
69 return eps;
70 }
71
72 [[nodiscard]] std::vector<std::string> get_irreps() const {
73 return irreps;
74 }
75
76 std::vector<int> get_localize_sets() const {
77 return localize_sets;
78 }
79
81 return occ;
82 }
83
84 /// setters will always invalidate all other member variables
85 MolecularOrbitals& set_mos(const std::vector<Function<T,NDIM> >& mo_new) {
87 mo=mo_new;
88 return *this;
89 }
90
91 /// updates will keep other member variables
92 MolecularOrbitals& update_mos(const std::vector<Function<T,NDIM> >& mo_new) {
93 mo=mo_new;
94 return *this;
95 }
96
98 occ=occ_new;
99 return *this;
100 }
101
102 /// updates will keep other member variables
103 MolecularOrbitals& update_localize_set(const std::vector<int>& set) {
104 localize_sets=set;
105 return *this;
106 }
107
108 /// updates will keep other member variables
110 const Tensor<double>& eps_new) {
111 mo=mo_new;
112 eps=copy(eps_new);
113 return *this;
114 }
115
116 MolecularOrbitals& recompute_irreps(const std::string pointgroup,
117 const Function<typename Tensor<T>::scalar_type,NDIM>& metric);
118
119 /// group orbitals into sets of similar orbital energies for localization
120 MolecularOrbitals& recompute_localize_sets(const double bandwidth=1.5) {
122 std::size_t nmo = mo.size();
123 std::vector<int> set = std::vector<int>(static_cast<size_t>(nmo), 0);
124 for (size_t i = 1; i < nmo; ++i) {
125 set[i] = set[i - 1];
126 // Only the new/boys localizers can tolerate not separating out the core orbitals
127 if (eps(i) - eps(i - 1) > bandwidth || get_occ()(i) != 1.0) ++(set[i]);
128 }
130 return *this;
131 }
132
133 static std::vector<Slice> convert_set_to_slice(const std::vector<int>& localized_set) {
134 std::vector<Slice> blocks;
135 long ilo=0;
136 for (size_t i=1; i<localized_set.size(); ++i) {
137 if (not (localized_set[i]==localized_set[i-1])) {
138 blocks.push_back(Slice(ilo, i-1));
139 ilo=i;
140 }
141 }
142 // add final block
143 blocks.push_back(Slice(ilo,localized_set.size()-1));
144 return blocks;
145 }
146
147
149 occ=Tensor<double>(mo.size());
150 occ=1.0;
151 return *this;
152 }
153
161
163 mo.clear();
164 }
165
167 eps.clear();
168 }
169
171 irreps.clear();
172 }
173
175 occ.clear();
176 }
177
179 localize_sets.clear();
180 }
181
182 void pretty_print(std::string message, std::vector<std::string> flags=std::vector<std::string>()) const {
183 print(message);
184 if (flags.size()==0) flags.resize(mo.size());
185 std::vector<std::string> irreps=get_irreps();
186 if (irreps.size()==0) irreps=std::vector<std::string>(mo.size(),"unknown");
187 print("orbital # irrep energy occupation localize_set");
188 for (int i=mo.size()-1; i>=0; --i) {
189// double n=get_mos()[i].norm2();
190 constexpr std::size_t bufsize=1024;
191 char buf[bufsize];
192 snprintf(buf,bufsize,"%5d %10s %12.8f %6.2f %8d %15s", i, irreps[i].c_str(),get_eps()[i],
193 get_occ()[i],get_localize_sets()[i], flags[i].c_str());
194 cout << std::string(buf) <<endl;
195 }
196 }
197
198
199 void print_frozen_orbitals(const long freeze) const {
200
201 World& world=mo.front().world();
202 MolecularOrbitals<T, 3> dummy_mo(*this);
203 dummy_mo.recompute_localize_sets();
204 if (world.rank() == 0) {
205 auto flags=std::vector<std::string>(dummy_mo.get_mos().size(),"active");
206 for (int i=0; i<freeze; ++i) flags[i]="frozen";
207 dummy_mo.pretty_print("diagonal Fock matrix elements with core/valence separation for freezing",flags);
208 print("\nfreezing orbitals: ", freeze,"\n");
209 }
210
211 }
212
213
214 /// @param[in] cubefile_header header of the cube file, from molecule::cubefile_header()
215 void print_cubefiles(const std::string name, const std::vector<std::string> cubefile_header) const;
216
217 template <typename Archive>
218 void serialize (Archive& ar) {
219 std::size_t nmo=mo.size();
220 ar & nmo;
221 if (nmo!=mo.size()) mo.resize(nmo);
222 for (auto& m : mo) ar & m;
223 ar & eps & irreps & occ & localize_sets;
224 if (ar.is_input_archive) {
225 if (irreps.size()==0) irreps=std::vector<std::string>(nmo,"unknown");
226 if (localize_sets.size()==0) localize_sets=std::vector<int>(nmo,0);
227 if (occ.size()==0) occ=Tensor<double>(nmo);
228 if (eps.size()==0) eps=Tensor<double>(nmo);
229 }
230 }
231
232 friend bool similar(const MolecularOrbitals& mo1, const MolecularOrbitals& mo2, const double thresh=1.e-6) {
233
234 if (mo1.mo.size()!=mo2.mo.size()) return false;
235 if (mo1.mo.size()==0) return true;
236
237 World& world=mo1.mo.front().world();
238 bool similar=((mo1.eps-mo2.eps).normf()<thresh);
239 similar=similar and (norm2(world,mo1.mo-mo2.mo)<thresh);
240 similar=similar and (mo1.irreps==mo2.irreps);
242 return similar;
243 }
244
245 void write_to(std::vector<Function<T,NDIM> >& mo_out, Tensor<double>& eps_out,
246 std::vector<std::string>& irrep_out, Tensor<double>& occ_out, std::vector<int>& set_out) const;
247
248 /// reads amo and bmo from the restartdata file
249
250 /// @return amo and bmo
251 std::pair<MolecularOrbitals<T,NDIM>, MolecularOrbitals<T,NDIM> >
252 static read_restartdata(World& world, const std::string filename, const Molecule& molecule,
253 const std::size_t nmo_alpha, const std::size_t nmo_beta) {
254 /*
255 * save and load from SCF.cc now contain the following extra variables
256 * To read Molecular Orbitals properly we need to take care of this data
257 * even if it is unused.
258 *
259 *
260 unsigned int version;
261 double L;
262 int k;
263 Molecule molecule;
264 std::string xc;
265
266*/
267 bool spinrestricted = false;
268 double L=0;
269 int k1=0; // Ignored for restarting, used in response only
270 unsigned int version = 4; // UPDATE THIS IF YOU CHANGE ANYTHING
271 // unsigned int archive_version;
272 double current_energy=0, converged_to_thresh=0;
273 std::string xc, localize_method;
274
275
277 ar & version;
278 ar & current_energy & spinrestricted;
279 ar & L& k1& molecule& xc & localize_method & converged_to_thresh;
280
282 amo.load_mos(ar, molecule, nmo_alpha);
283 bool have_beta=(not spinrestricted) and (nmo_beta>0);
284 if (have_beta) {
285 bmo.load_mos(ar,molecule,nmo_beta);
286 }
287 return std::make_pair(amo,bmo);
288 }
289
290 /// reads amo and bmo from the restartdata file
291
292 /// @return amo and bmo
293 void static save_restartdata(World& world, const std::string filename, const Molecule& molecule,
295 // TODO there is a good chance that this needs to be modified if it is intended to be read by SCF save/load
296 bool spinrestricted = false;
297 unsigned int version=4;
298 double current_energy=0.0, converged_to_thresh=1.e10;
299 double L=0;
300 std::string xc, localize_method;
301 int k1=0; // Ignored for restarting, used in response only
303 ar & version;
304 ar & current_energy & spinrestricted;
305 ar & L& k1& molecule& xc & localize_method & converged_to_thresh;
306
307 amo.save_mos(ar,molecule);
308 bmo.save_mos(ar,molecule);
309 }
310
311 /// legacy code
312 void load_mos(archive::ParallelInputArchive<>& ar, const Molecule& molecule, const std::size_t nmo_from_input) {
313
314 unsigned int nmo = 0;
315
316 ar & nmo;
317 MADNESS_ASSERT(nmo >= nmo_from_input);
318 ar & eps & occ & localize_sets;
319 mo.resize(nmo);
320 for (unsigned int i = 0; i < mo.size(); ++i)
321 ar & mo[i];
322 unsigned int n_core = molecule.n_core_orb_all();
323 if (nmo > nmo_from_input) {
324 localize_sets = vector<int>(localize_sets.begin() + n_core,
325 localize_sets.begin() + n_core + nmo_from_input);
326 mo = std::vector<Function<T,NDIM> >(mo.begin() + n_core,
327 mo.begin() + n_core + nmo_from_input);
328 eps = copy(eps(Slice(n_core, n_core + nmo_from_input - 1)));
329 occ = copy(occ(Slice(n_core, n_core + nmo_from_input - 1)));
330 }
331 }
332
333 /// legacy code
335
336 unsigned int nmo=mo.size();
337
338 ar & nmo;
339 ar & eps & occ & localize_sets;
340 for (unsigned int i = 0; i < mo.size(); ++i)
341 ar & mo[i];
342 //unsigned int n_core = molecule.n_core_orb_all();
343 }
344
345 void post_process_mos(World& world, const double thresh, const int k);
346
347
348 /// save MOs in the AO projection for geometry restart
349
350 /// compute the aos in MRA projection as:
351 static void save_restartaodata(World& world, const Molecule& molecule,
353 const AtomicBasisSet& aobasis);
354
355 /// uses AO-projection as a restart guess
356
357 /// @return amo and bmo
358 std::pair<MolecularOrbitals<T,NDIM>, MolecularOrbitals<T,NDIM> >
359 static read_restartaodata(World& world,
360 const Molecule& molecule, const bool have_beta);
361
362 void project_ao(World& world, const Tensor<T>& Saomo, const std::vector<Function<double,3> >& aos);
363
364 std::vector<Vector<typename Tensor<T>::scalar_type,3>> compute_center(
365 const Function<typename Tensor<T>::scalar_type,NDIM> metric2=Function<typename Tensor<T>::scalar_type,NDIM>()) const;
366
367private:
368 std::vector<Function<T,NDIM> > mo;
370 std::vector<std::string> irreps;
372 std::vector<int> localize_sets;
373
374};
375
376} /* namespace madness */
377
378#endif /* SRC_APPS_CHEM_MOLECULARORBITALS_H_ */
Contracted Gaussian basis.
Definition apps/periodic_old/molecularbasis.h:424
Contracted Gaussian basis.
Definition madness/chem/molecularbasis.h:465
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
A multiresolution adaptive numerical function.
Definition mra.h:122
Definition MolecularOrbitals.h:24
void invalidate_eps()
Definition MolecularOrbitals.h:166
MolecularOrbitals get_subset(const int iset) const
Definition MolecularOrbitals.h:45
static void save_restartaodata(World &world, const Molecule &molecule, const MolecularOrbitals< T, NDIM > &amo, const MolecularOrbitals< T, NDIM > &bmo, const AtomicBasisSet &aobasis)
save MOs in the AO projection for geometry restart
Definition MolecularOrbitals.cc:83
std::vector< std::string > irreps
Definition MolecularOrbitals.h:370
void print_cubefiles(const std::string name, const std::vector< std::string > cubefile_header) const
Definition MolecularOrbitals.cc:66
void write_to(std::vector< Function< T, NDIM > > &mo_out, Tensor< double > &eps_out, std::vector< std::string > &irrep_out, Tensor< double > &occ_out, std::vector< int > &set_out) const
Definition MolecularOrbitals.cc:34
static std::pair< MolecularOrbitals< T, NDIM >, MolecularOrbitals< T, NDIM > > read_restartdata(World &world, const std::string filename, const Molecule &molecule, const std::size_t nmo_alpha, const std::size_t nmo_beta)
reads amo and bmo from the restartdata file
Definition MolecularOrbitals.h:252
MolecularOrbitals & update_mos(const std::vector< Function< T, NDIM > > &mo_new)
updates will keep other member variables
Definition MolecularOrbitals.h:92
MolecularOrbitals & recompute_localize_sets(const double bandwidth=1.5)
group orbitals into sets of similar orbital energies for localization
Definition MolecularOrbitals.h:120
Tensor< double > get_eps() const
Definition MolecularOrbitals.h:68
std::size_t size() const
Definition MolecularOrbitals.h:60
void save_mos(archive::ParallelOutputArchive<> &ar, const Molecule &molecule) const
legacy code
Definition MolecularOrbitals.h:334
void serialize(Archive &ar)
Definition MolecularOrbitals.h:218
friend bool similar(const MolecularOrbitals &mo1, const MolecularOrbitals &mo2, const double thresh=1.e-6)
Definition MolecularOrbitals.h:232
static std::pair< MolecularOrbitals< T, NDIM >, MolecularOrbitals< T, NDIM > > read_restartaodata(World &world, const Molecule &molecule, const bool have_beta)
uses AO-projection as a restart guess
Definition MolecularOrbitals.cc:104
void invalidate_occ()
Definition MolecularOrbitals.h:174
Tensor< double > eps
Definition MolecularOrbitals.h:369
MolecularOrbitals(const std::vector< Function< T, NDIM > > &mo)
Definition MolecularOrbitals.h:31
static std::vector< Slice > convert_set_to_slice(const std::vector< int > &localized_set)
Definition MolecularOrbitals.h:133
void invalidate_mos()
Definition MolecularOrbitals.h:162
MolecularOrbitals & set_all_orbitals_occupied()
Definition MolecularOrbitals.h:148
void invalidate_irreps()
Definition MolecularOrbitals.h:170
void invalidate_localize_sets()
Definition MolecularOrbitals.h:178
MolecularOrbitals & update_occ(const Tensor< double > &occ_new)
Definition MolecularOrbitals.h:97
MolecularOrbitals & update_localize_set(const std::vector< int > &set)
updates will keep other member variables
Definition MolecularOrbitals.h:103
MolecularOrbitals & update_mos_and_eps(const std::vector< Function< T, NDIM > > &mo_new, const Tensor< double > &eps_new)
updates will keep other member variables
Definition MolecularOrbitals.h:109
MolecularOrbitals(const std::vector< Function< T, NDIM > > &mo, const Tensor< double > &eps)
Definition MolecularOrbitals.h:35
void pretty_print(std::string message, std::vector< std::string > flags=std::vector< std::string >()) const
Definition MolecularOrbitals.h:182
MolecularOrbitals(const std::vector< Function< T, NDIM > > &mo, const Tensor< double > &eps, const std::vector< std::string > &irrep, const Tensor< double > &occ, const std::vector< int > &set)
Definition MolecularOrbitals.h:39
void project_ao(World &world, const Tensor< T > &Saomo, const std::vector< Function< double, 3 > > &aos)
Definition MolecularOrbitals.cc:141
Tensor< double > get_occ() const
Definition MolecularOrbitals.h:80
std::vector< int > get_localize_sets() const
Definition MolecularOrbitals.h:76
std::vector< Function< T, NDIM > > mo
Definition MolecularOrbitals.h:368
std::vector< Function< T, NDIM > > get_mos() const
Definition MolecularOrbitals.h:64
MolecularOrbitals(const MolecularOrbitals< T, NDIM > &other)=default
std::vector< int > localize_sets
Definition MolecularOrbitals.h:372
Tensor< double > occ
Definition MolecularOrbitals.h:371
std::vector< std::string > get_irreps() const
Definition MolecularOrbitals.h:72
static void save_restartdata(World &world, const std::string filename, const Molecule &molecule, const MolecularOrbitals< T, NDIM > &amo, const MolecularOrbitals< T, NDIM > &bmo)
reads amo and bmo from the restartdata file
Definition MolecularOrbitals.h:293
void invalidate_all()
Definition MolecularOrbitals.h:154
void print_frozen_orbitals(const long freeze) const
Definition MolecularOrbitals.h:199
void post_process_mos(World &world, const double thresh, const int k)
Definition MolecularOrbitals.cc:53
void load_mos(archive::ParallelInputArchive<> &ar, const Molecule &molecule, const std::size_t nmo_from_input)
legacy code
Definition MolecularOrbitals.h:312
std::vector< Vector< typename Tensor< T >::scalar_type, 3 > > compute_center(const Function< typename Tensor< T >::scalar_type, NDIM > metric2=Function< typename Tensor< T >::scalar_type, NDIM >()) const
Definition MolecularOrbitals.cc:155
MolecularOrbitals & set_mos(const std::vector< Function< T, NDIM > > &mo_new)
setters will always invalidate all other member variables
Definition MolecularOrbitals.h:85
MolecularOrbitals & recompute_irreps(const std::string pointgroup, const Function< typename Tensor< T >::scalar_type, NDIM > &metric)
Definition MolecularOrbitals.cc:24
Definition molecule.h:124
A slice defines a sub-range or patch of a dimension.
Definition slice.h:103
A tensor is a multidimension array.
Definition tensor.h:317
TensorTypeData< T >::scalar_type scalar_type
C++ typename of the real type associated with a complex type.
Definition tensor.h:409
void clear()
Frees all memory and resests to state of default constructor.
Definition tensor.h:1884
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
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
Objects that implement their own parallel archive interface should derive from this class.
Definition parallel_archive.h:58
const std::size_t bufsize
Definition derivatives.cc:16
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:182
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
static const char * filename
Definition legendre.cc:96
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 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
std::vector< std::string > cubefile_header(std::string filename="input", const bool &no_orient=false)
Definition molecule.cc:69
std::string name(const FuncType &type, const int ex=-1)
Definition ccpairfunction.h:28
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition mra.h:2002
XCfunctional xc
Definition newsolver_lda.cc:53
Implements ParallelInputArchive and ParallelOutputArchive for parallel serialization of data.
static const double m
Definition relops.cc:9
static const double L
Definition rk.cc:46
static const double thresh
Definition rk.cc:45
static const long k
Definition rk.cc:44
Defines and implements most of Tensor.
static const std::size_t NDIM
Definition testpdiff.cc:42
double k1
Definition testperiodic.cc:67
static Molecule molecule
Definition testperiodicdft.cc:38
static AtomicBasisSet aobasis
Definition testperiodicdft.cc:39