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 
17 namespace madness {
18 
19 template<typename T, std::size_t NDIM> class Function;
20 //class World;
21 class AtomicBasisSet;
22 
23 template<typename T, std::size_t NDIM>
25 public:
26 
28 
29  MolecularOrbitals() = default;
30 
31  MolecularOrbitals(const std::vector<Function<T,NDIM> >& mo)
32  : mo(mo), eps(), irreps(), occ(), localize_sets() {
33  }
34 
35  MolecularOrbitals(const std::vector<Function<T,NDIM> >& mo, const Tensor<double>& eps)
36  : mo(mo), eps(eps), irreps(), occ(), localize_sets() {
37  }
38 
39  MolecularOrbitals(const std::vector<Function<T,NDIM> >& mo, const Tensor<double>& eps,
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 {
46  auto slices = convert_set_to_slice(localize_sets);
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  }
129  update_localize_set(set);
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 
154  void invalidate_all() {
155  invalidate_mos();
156  invalidate_eps();
158  invalidate_occ();
160  }
161 
162  void invalidate_mos() {
163  mo.clear();
164  }
165 
166  void invalidate_eps() {
167  eps.clear();
168  }
169 
171  irreps.clear();
172  }
173 
174  void invalidate_occ() {
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 
281  MolecularOrbitals<T,NDIM> amo, bmo;
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,
294  const MolecularOrbitals<T,NDIM>& amo, const MolecularOrbitals<T,NDIM>& bmo) {
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 
367 private:
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
MolecularOrbitals & set_mos(const std::vector< Function< T, NDIM > > &mo_new)
setters will always invalidate all other member variables
Definition: MolecularOrbitals.h:85
std::vector< Function< T, NDIM > > get_mos() const
Definition: MolecularOrbitals.h:64
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
MolecularOrbitals & update_localize_set(const std::vector< int > &set)
updates will keep other member variables
Definition: MolecularOrbitals.h:103
MolecularOrbitals & set_all_orbitals_occupied()
Definition: MolecularOrbitals.h:148
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
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
Tensor< double > get_occ() const
Definition: MolecularOrbitals.h:80
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
std::vector< std::string > get_irreps() const
Definition: MolecularOrbitals.h:72
static std::vector< Slice > convert_set_to_slice(const std::vector< int > &localized_set)
Definition: MolecularOrbitals.h:133
MolecularOrbitals & update_occ(const Tensor< double > &occ_new)
Definition: MolecularOrbitals.h:97
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
void invalidate_mos()
Definition: MolecularOrbitals.h:162
std::vector< int > get_localize_sets() const
Definition: MolecularOrbitals.h:76
void invalidate_irreps()
Definition: MolecularOrbitals.h:170
void invalidate_localize_sets()
Definition: MolecularOrbitals.h:178
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
std::vector< Function< T, NDIM > > mo
Definition: MolecularOrbitals.h:368
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 MolecularOrbitals< T, NDIM > &other)=default
std::vector< int > localize_sets
Definition: MolecularOrbitals.h:372
Tensor< double > occ
Definition: MolecularOrbitals.h:371
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 & 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
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
const double m
Definition: gfit.cc:199
#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:190
#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
const char * version()
Get the MADNESS version number.
Definition: info.cc:17
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
static const char * filename
Definition: legendre.cc:96
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
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
XCfunctional xc
Definition: newsolver_lda.cc:53
Implements ParallelInputArchive and ParallelOutputArchive for parallel serialization of data.
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