MPQC  3.0.0-alpha
findisp.h
1 //
2 // findisp.h
3 //
4 // Copyright (C) 1997 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifndef _chemistry_molecule_findisp_h
29 #define _chemistry_molecule_findisp_h
30 
31 #include <iostream>
32 #include <map>
33 #include <vector>
34 #include <algorithm>
35 
36 #include <util/misc/scexception.h>
37 #include <chemistry/molecule/deriv.h>
38 #include <chemistry/molecule/energy.h>
39 
40 namespace sc {
41 
44 
47  template <typename Value>
48  class Displacements {
49  public:
50  Displacements() {}
51  ~Displacements() {}
52 
53  void save(StateOut& s) const {
54  s.put(disps_);
55  }
56  void restore(StateIn& s) {
57  s.get(disps_);
58  }
59 
60  typedef int Coord;
61  typedef double DisplacementSize;
62  typedef std::pair<Coord,DisplacementSize> KeyElem;
63  typedef std::vector<KeyElem> Key;
64  typedef typename std::map< Key, Value>::iterator iter;
65  typedef typename std::map< Key, Value>::const_iterator citer;
66 
67  citer begin() const { return disps_.begin(); }
68  citer end() const { return disps_.end(); }
69  std::size_t size() const { return disps_.size(); }
70 
71  std::pair<Key,Value> find(Coord c1, DisplacementSize d1) const {
72  Key key;
73  key.push_back(std::make_pair(c1,d1));
74  return find(key);
75  }
76  std::pair<Key,Value> find(Coord c1, DisplacementSize d1,
77  Coord c2, DisplacementSize d2) const {
78  Key key;
79  key.push_back(std::make_pair(c1,d1));
80  key.push_back(std::make_pair(c2,d2));
81  std::stable_sort(key.begin(), key.end());
82  return find(key);
83  }
85  bool has(const Key& key) const {
86  citer v = disps_.find(key);
87  return v != disps_.end();
88  }
90  std::pair<Key,Value> find(const Key& key) const {
91  citer v = disps_.find(key);
92  if (v == disps_.end())
93  throw sc::ProgrammingError("Displacement::find -- displacement not found",__FILE__,__LINE__);
94  return *v;
95  }
97  citer find(const Value& value) const {
98  for(citer v = disps_.begin(); v != disps_.end(); ++v) {
99  if ((*v).second == value) return v;
100  }
101  throw sc::ProgrammingError("Displacement::find -- value not found",__FILE__,__LINE__);
102  }
103 
104  void push(Coord c1, DisplacementSize d1,
105  const Value& v) {
106  Key key;
107  key.push_back(std::make_pair(c1,d1));
108  return push(key,v);
109  }
110  void push(Coord c1, DisplacementSize d1,
111  Coord c2, DisplacementSize d2,
112  const Value& v) {
113  Key key;
114  key.push_back(std::make_pair(c1,d1));
115  key.push_back(std::make_pair(c2,d2));
116  std::stable_sort(key.begin(), key.end());
117  return push(key,v);
118  }
120  void push(const Key& key, const Value& v) {
121  iter i = disps_.find(key);
122  if (i != disps_.end())
123  throw sc::ProgrammingError("Displacement::push -- displacement already exists",__FILE__,__LINE__);
124  disps_[key] = v;
125  }
126 
127  private:
133  std::map< Key, Value> disps_;
134  };
135 
136 #if 0
137 
141  template <unsigned int TargetOrder, typename Function>
142  class FinDispDerivative : virtual public SavableState {
143  public:
149  FinDispDerivative(const Ref<Function>& function);
150  FinDispDerivative(StateIn&);
151  ~FinDispDerivative();
152  void save_data_state(StateOut&);
153 
154  typedef typename Function::Value FunctionValue;
155  typedef typename Function::Parameter FunctionParameter;
156  static const unsigned int TensorTraits<FunctionParameter>::rank param_rank;
157  const ParameterDerivative<FunctionValue,IntegerProduct<TargetOrder,param_rank>::value>::result Result;
158 
159  const Result& result();
160  private:
161  static ClassDesc class_desc_;
162  Ref<Function> function_;
163  Result result_;
164 
165  void compute();
166  };
167 #endif
168 
170  struct EGH {
171  public:
172  EGH();
173  EGH(double e, const RefSCVector& g, const RefSymmSCMatrix& h);
174  ~EGH();
175 
176  double energy() const { return energy_; }
177  const RefSCVector& gradient() const { return gradient_; }
178  const RefSymmSCMatrix& hessian() const { return hessian_; }
179 
180  private:
181  double energy_;
182  RefSCVector gradient_;
183  RefSymmSCMatrix hessian_;
184  };
186  template <> void FromStateIn<EGH>(EGH& v, StateIn& s, int& count);
187  template <> void ToStateOut<EGH>(const EGH& v, StateOut& so, int& count);
189 
190 
195 
197  class Params : virtual public SavableState {
198  public:
199  Params();
200  Params(const Ref<KeyVal>& kv);
201  Params(StateIn&);
202  ~Params();
203  void save_data_state(StateOut&);
204 
205  const Ref<PointGroup>& disp_pg() const { return disp_pg_; }
206  double disp_size() const { return disp_; }
207  bool only_totally_symmetric() const { return only_totally_symmetric_; }
208  bool user_provided_eliminate_quadratic_terms() const { return user_provided_eliminate_quadratic_terms_; }
209  bool eliminate_quadratic_terms() const { return eliminate_quadratic_terms_; }
210  bool do_null_displacement() const { return do_null_displacement_; }
211  int debug() const { return debug_; }
212  bool checkpoint() const { return checkpoint_; }
213  const std::string& checkpoint_file() const { return checkpoint_file_; }
214  bool restart() const { return restart_; }
215  const std::string& restart_file() const { return restart_file_; }
216  bool use_energies() const { return use_energies_; }
217  double gradient_accuracy() const { return gradient_accuracy_; }
218  double energy_accuracy() const { return energy_accuracy_; }
219  int nirrep() const { return disp_pg_->char_table().nirrep(); }
220 
221  // if called will no need to (re)compute default
222  void set_eliminate_quadratic_terms(bool e) {
223  user_provided_eliminate_quadratic_terms_ = true;
224  eliminate_quadratic_terms_ = e;
225  }
226  void set_disp_size(double s);
227  void set_disp_pg(const Ref<PointGroup>& pg) { disp_pg_ = pg; }
228  void set_restart(bool r = true) { restart_ = r; }
229  void set_checkpoint(bool c = true) { checkpoint_ = c; }
230  void set_desired_accuracy(double acc);
231 
232  private:
233  static ClassDesc class_desc_;
234 
235  // In case molecule must be given in lower symmetry, its actual
236  // symmetry and the symmetry used to compute displacements is this
237  Ref<PointGroup> disp_pg_;
238  // the cartesian displacement size in bohr
239  double disp_;
240  // only do the totally symmetric displacements (makes sense if the Hessian is to be used in geometry optimization)
241  bool only_totally_symmetric_;
242  // eliminate the quadratic terms in the hessian by doing an extra displacement for
243  // each of the totally-symmetric coordinates
244  bool eliminate_quadratic_terms_;
245  bool user_provided_eliminate_quadratic_terms_; // whether user provided the value eliminate_quadratic_terms
246  // if true, will skip the MolecularEnergy-dependent default override
247  // use the gradient at the initial geometry to remove first order terms
248  // (important if not at equilibrium geometry)
249  bool do_null_displacement_;
250  // print flag
251  int debug_;
252  // whether or not to checkpoint
253  bool checkpoint_;
254  // the name of the checkpoint file
255  std::string checkpoint_file_;
256  // whether or not to attempt a restart
257  bool restart_;
258  // the name of the restart file
259  std::string restart_file_;
260  // force computation from energies
261  bool use_energies_;
262  // the accuracy for energy calculations
263  double energy_accuracy_;
264  // the accuracy for gradient calculations
265  double gradient_accuracy_;
266  };
267 
268 
269  // Implements FinDispMolecularHessian
270  class Impl : virtual public SavableState {
271  public:
272  typedef Displacements<EGH>::Key Displacement;
273 
274  Impl(const Ref<MolecularEnergy>& e,
275  const Ref<Params>& params);
276  Impl(StateIn& s);
277  ~Impl();
278  void save_data_state(StateOut& s);
279 
280  void checkpoint_displacements(StateOut&);
281  void restore_displacements(StateIn&);
282 
283  const Ref<Params>& params() const { return params_; }
284  const Ref<MolecularEnergy>& mole() const { return mole_; }
285  void set_mole(const Ref<MolecularEnergy>& mole) { mole_ = mole; }
286  virtual int ndisplace() const =0;
287  // returns a matrix whose columns are SALC displacements
288  RefSCMatrix displacements(int irrep) const;
289  void displace(const Displacement& d);
290  void original_geometry();
291 
292  virtual void init();
293  virtual void restart();
294 
298 
299  // transforms hessian in symmetrized coordinates to cartesian coordinates
300  void do_hess_for_irrep(int irrep,
301  const RefSymmSCMatrix &dhessian,
302  const RefSymmSCMatrix &xhessian);
303 
304  protected:
305  Ref<Params> params_;
306 
307  Ref<MolecularEnergy> mole_;
308  // The molecule's original point group for restoration at the end.
309  Ref<PointGroup> original_point_group_;
310  // The molecule's original geometry for restoration at the end and
311  //computing displacements.
312  RefSCVector original_geometry_;
313  // a basis for the symmetrized cartesian coordinates (rows)
314  RefSCMatrix symbasis_;
315 
316  // values, gradients, and hessians at various (displaced) geometries
317  Displacements<EGH> values_;
318 
320  virtual RefSymmSCMatrix compute_hessian() =0;
322  virtual void compute_mole(const Displacement& d) =0;
323 
324  Ref<SCMatrixKit> matrixkit() const { return mole_->matrixkit(); }
325  RefSCDimension d3natom() const { return mole_->moldim(); }
327  unsigned int coor_to_irrep(unsigned int symm_coord) const;
328 
329  static ClassDesc class_desc_;
330  };
331 
332  // Implementation of finite-difference hessians from gradients
333  class GradientsImpl : public Impl {
334  public:
335  GradientsImpl(const Ref<MolecularEnergy>& e,
336  const Ref<Params>& params);
337  GradientsImpl(StateIn& s);
338  ~GradientsImpl();
339  void save_data_state(StateOut& s);
340 
341  private:
342 
343  // will throw if mole lacks gradients capability; do nothing if mole is null
344  static void validate_mole(const Ref<MolecularEnergy>& e);
345  int ndisplace() const;
346  void set_gradient(const Displacement& d, double energy, const RefSCVector &grad);
347  RefSymmSCMatrix compute_hessian();
348  void compute_mole(const Displacement& d);
349 
350  static ClassDesc class_desc_;
351  };
352 
353  // Implementation of finite-difference hessians from energies
354  class EnergiesImpl : public Impl {
355  public:
356  // will throw if mole lacks energies capability; do nothing if mole is null
357  EnergiesImpl(const Ref<MolecularEnergy>& e,
358  const Ref<Params>& params);
359  EnergiesImpl(StateIn& s);
360  ~EnergiesImpl();
361  void save_data_state(StateOut& s);
362 
363  private:
364 
365  static void validate_mole(const Ref<MolecularEnergy>& e);
366  int ndisplace() const;
367  RefSymmSCMatrix compute_hessian();
368  void compute_mole(const Displacement& d);
369  const Displacements<EGH>& values() const { return values_; }
370 
371  struct Eij {
372  Eij(int i, int j, EnergiesImpl& eimpl) :
373  i_(i), j_(j), eimpl_(eimpl)
374  {
375  }
376 
377  // return energy computed at geometry displaced by di (in units of disp_)
378  // along i and dj along j
379  double operator()(int di, int dj) {
380  Displacement disp;
381  if (i_ != j_) {
382  if (di != 0) disp.push_back(std::make_pair(i_,double(di)));
383  if (dj != 0) disp.push_back(std::make_pair(j_,double(dj)));
384  }
385  else { // i_ == j_
386  if (di+dj != 0) disp.push_back(std::make_pair(i_,double(di+dj)));
387  }
388  eimpl_.compute_mole(disp);
389  const double energy = eimpl_.values().find(disp).second.energy();
390  return energy;
391  }
392 
393  int i_, j_;
394  EnergiesImpl& eimpl_;
395  };
396 
397  static ClassDesc class_desc_;
398  };
399 
400  private:
401  Ref<Impl> pimpl_; //<< initliazed lazily
402  Ref<MolecularEnergy> mole_init_; //< pimpl_ is initalized lazily, so this is used to hold the MolecularEnergy object used to initialize the pimpl_
403  Ref<Params> params_;
404 
405  void override_default_params(); // override some defaults based on the properties of MolecularEnergy
406 
407  protected:
408 
410  void init_pimpl(const Ref<MolecularEnergy>& e);
411  void restart();
412 
413  public:
485  void save_data_state(StateOut&);
486 
490 
492  MolecularEnergy* energy() const { return pimpl_ ? pimpl_->mole().pointer() : 0; }
493 
494  const Ref<Params>& params() const { return params_; }
495 
496  void set_desired_accuracy(double acc);
497 };
498 
503 //
504 // public:
505 //
506 // class Displacements {
507 // public:
508 //
509 // private:
510 // /** displacements are represented as a key -> index map
511 // * keys are formatted as comma-separated list of integer,double pairs
512 // * each pair contains the index of SALC (0-based) + displacement size
513 // * SALC indices are nondecreasing
514 // */
515 // std::map< std::list< std::pair<int,double> >, > disps;
516 // };
517 
518  protected:
519  Ref<MolecularEnergy> mole_;
520  // In case molecule must be given in lower symmetry, its actual
521  // symmetry and the symmetry used to compute displacements is this
522  Ref<PointGroup> displacement_point_group_;
523  // The molecule's original geometry for restoration at the end and
524  //computing displacements.
525  RefSCVector original_geometry_;
526  // the cartesian displacement size in bohr
527  double disp_;
528  // the accuracy of the energy computations
529  double energy_accuracy_;
530  // whether or not to attempt a restart
531  int restart_;
532  // the name of the restart file
533  std::string restart_file_;
534  // whether or not to checkpoint
535  int checkpoint_;
536  // the name of the checkpoint file
537  std::string checkpoint_file_;
538  // 2-pt formula for gradient is accurate to O(h^2) (contaminated by 3rd derivatives)
539  // make it O(h^4), i.e. eliminate the h^2 terms, by using a 4-pt formula
540  int eliminate_quadratic_terms_;
541  // print flag
542  int debug_;
543  // a basis for the symmetrized cartesian coordinates
544  RefSCMatrix symbasis_;
545  // the energies at each of the completed displacements
546  std::vector<double> energies_;
547 
548  // given displacement # disp returns internal coordinate index and displacement size (in units of disp_)
549  void get_disp(int disp, int &index, double &dispsize);
550 // void do_grad_for_irrep(int irrep,
551 // const RefSymmSCMatrix &dhessian,
552 // const RefSymmSCMatrix &xhessian);
553  void init();
554  void restart();
555 
558  RefSCVector compute_gradient();
559  int ndisplace() const;
560  int ndisplacements_done() const { return energies_.size(); }
561  RefSCMatrix displacements(int irrep) const;
562  void displace(int disp);
563  void original_geometry();
564  //void set_gradient(int disp, const RefSCVector &grad);
565  void checkpoint_displacements(StateOut&);
566  void restore_displacements(StateIn&);
567 
568  public:
619  void save_data_state(StateOut&);
620 
621 
625 
627  void set_checkpoint(int c) { checkpoint_ = c; }
629  int checkpoint() const { return checkpoint_; }
630 
632  MolecularEnergy* energy() const;
633 
634  Ref<SCMatrixKit> matrixkit() const { return mole_->matrixkit(); }
635  RefSCDimension d3natom() const { return mole_->moldim(); }
636 
637  void set_desired_accuracy(double acc);
638 
639  void set_eliminate_quadratic_terms(bool e) { eliminate_quadratic_terms_ = e; }
640  void set_disp_size(double s);
641 };
642 
644 // end of addtogroup ChemistryMolecule
645 
646 }
647 
648 #endif
649 
650 // Local Variables:
651 // mode: c++
652 // c-file-style: "CLJ"
653 // End:
sc::FinDispMolecularGradient::save_data_state
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
sc::FinDispMolecularGradient
Computes the molecular gradient by finite differences of energies.
Definition: findisp.h:502
sc::RefSymmSCMatrix
The RefSymmSCMatrix class is a smart pointer to an SCSymmSCMatrix specialization.
Definition: matrix.h:265
sc::FinDispMolecularGradient::checkpoint
int checkpoint() const
Return the current value of the checkpoint option.
Definition: findisp.h:629
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition: matrix.h:135
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::FinDispMolecularHessian::set_desired_accuracy
void set_desired_accuracy(double acc)
Sets the desired accuracy.
sc::Ref::pointer
T * pointer() const
Returns a pointer the reference counted object.
Definition: ref.h:413
sc::MolecularGradient
MolecularGradient is an abstract class that computes a molecule's first derivatives of the energy wit...
Definition: deriv.h:207
sc::MolecularHessian
MolecularHessian is an abstract class that computes a molecule's second derivatives of the energy wit...
Definition: deriv.h:46
sc::Displacements::has
bool has(const Key &key) const
assumes that key is sorted
Definition: findisp.h:85
sc::FinDispMolecularHessian::init_pimpl
void init_pimpl(const Ref< MolecularEnergy > &e)
initializes pimpl_, it should not be called until e is fully initalized, hence use this lazily
sc::StateIn
Definition: statein.h:79
sc::RefSCDimension
The RefSCDimension class is a smart pointer to an SCDimension specialization.
Definition: dim.h:152
sc::ClassDesc
This class is used to contain information about classes.
Definition: class.h:147
sc::FinDispMolecularHessian::set_energy
void set_energy(const Ref< MolecularEnergy > &energy)
Some MolecularHessian specializations require a molecular energy object.
sc::EGH
energy + gradient + hessian
Definition: findisp.h:170
sc::MolecularEnergy::energy
virtual double energy()
A wrapper around value();.
sc::MolecularEnergy
The MolecularEnergy abstract class inherits from the Function class.
Definition: energy.h:50
sc::FinDispMolecularGradient::set_checkpoint
void set_checkpoint(int c)
Set checkpoint option.
Definition: findisp.h:627
sc::FinDispMolecularGradient::set_desired_accuracy
void set_desired_accuracy(double acc)
Sets the desired accuracy.
sc::FinDispMolecularHessian::save_data_state
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
sc::RefSCVector
The RefSCVector class is a smart pointer to an SCVector specialization.
Definition: matrix.h:55
sc::StateOut
Definition: stateout.h:71
sc::FinDispMolecularHessian
Computes the molecular hessian by finite displacements of gradients (or, if not available,...
Definition: findisp.h:194
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition: scexception.h:92
sc::FinDispMolecularGradient::cartesian_gradient
RefSCVector cartesian_gradient()
This returns the cartesian gradient.
sc::FinDispMolecularGradient::set_energy
void set_energy(const Ref< MolecularEnergy > &energy)
Some MolecularGradient specializations require a molecular energy object.
sc::StateOut::put
virtual int put(const ClassDesc *)
Write out information about the given ClassDesc.
sc::SavableState
Base class for objects that can save/restore state.
Definition: state.h:45
sc::Displacements::find
citer find(const Value &value) const
finds the first instance of value (implemented as dumb O(N) search)
Definition: findisp.h:97
sc::StateIn::get
virtual int get(const ClassDesc **)
This restores ClassDesc's.
sc::Displacements::push
void push(const Key &key, const Value &v)
assumes that key is sorted
Definition: findisp.h:120
sc::FinDispMolecularHessian::energy
MolecularEnergy * energy() const
This returns a MolecularEnergy object, if used by this specialization.
Definition: findisp.h:492
sc::FinDispMolecularHessian::cartesian_hessian
RefSymmSCMatrix cartesian_hessian()
This returns the cartesian hessian.
sc::Displacements::find
std::pair< Key, Value > find(const Key &key) const
assumes that key is sorted
Definition: findisp.h:90
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::Displacements
Maps displacements in terms of symmetrized coordinates to property values.
Definition: findisp.h:48
sc::FinDispMolecularGradient::energy
MolecularEnergy * energy() const
This returns a MolecularEnergy object, if used by this specialization.

Generated at Sun Jan 26 2020 23:23:56 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.