28 #ifndef _chemistry_molecule_findisp_h
29 #define _chemistry_molecule_findisp_h
36 #include <util/misc/scexception.h>
37 #include <chemistry/molecule/deriv.h>
38 #include <chemistry/molecule/energy.h>
47 template <
typename Value>
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;
67 citer begin()
const {
return disps_.begin(); }
68 citer end()
const {
return disps_.end(); }
69 std::size_t size()
const {
return disps_.size(); }
71 std::pair<Key,Value> find(Coord c1, DisplacementSize d1)
const {
73 key.push_back(std::make_pair(c1,d1));
76 std::pair<Key,Value> find(Coord c1, DisplacementSize d1,
77 Coord c2, DisplacementSize d2)
const {
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());
85 bool has(
const Key& key)
const {
86 citer v = disps_.find(key);
87 return v != disps_.end();
90 std::pair<Key,Value>
find(
const Key& key)
const {
91 citer v = disps_.find(key);
92 if (v == disps_.end())
97 citer
find(
const Value& value)
const {
98 for(citer v = disps_.begin(); v != disps_.end(); ++v) {
99 if ((*v).second == value)
return v;
104 void push(Coord c1, DisplacementSize d1,
107 key.push_back(std::make_pair(c1,d1));
110 void push(Coord c1, DisplacementSize d1,
111 Coord c2, DisplacementSize d2,
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());
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__);
133 std::map< Key, Value> disps_;
141 template <
unsigned int TargetOrder,
typename Function>
142 class FinDispDerivative :
virtual public SavableState {
149 FinDispDerivative(
const Ref<Function>&
function);
150 FinDispDerivative(StateIn&);
151 ~FinDispDerivative();
152 void save_data_state(StateOut&);
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;
159 const Result& result();
161 static ClassDesc class_desc_;
162 Ref<Function> function_;
176 double energy()
const {
return energy_; }
177 const RefSCVector& gradient()
const {
return gradient_; }
186 template <>
void FromStateIn<EGH>(
EGH& v,
StateIn& s,
int& count);
187 template <>
void ToStateOut<EGH>(
const EGH& v,
StateOut& so,
int& count);
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(); }
222 void set_eliminate_quadratic_terms(
bool e) {
223 user_provided_eliminate_quadratic_terms_ =
true;
224 eliminate_quadratic_terms_ = e;
226 void set_disp_size(
double s);
228 void set_restart(
bool r =
true) { restart_ = r; }
229 void set_checkpoint(
bool c =
true) { checkpoint_ = c; }
241 bool only_totally_symmetric_;
244 bool eliminate_quadratic_terms_;
245 bool user_provided_eliminate_quadratic_terms_;
249 bool do_null_displacement_;
255 std::string checkpoint_file_;
259 std::string restart_file_;
263 double energy_accuracy_;
265 double gradient_accuracy_;
272 typedef Displacements<EGH>::Key Displacement;
280 void checkpoint_displacements(
StateOut&);
281 void restore_displacements(
StateIn&);
283 const Ref<Params>& params()
const {
return params_; }
286 virtual int ndisplace()
const =0;
289 void displace(
const Displacement& d);
290 void original_geometry();
293 virtual void restart();
300 void do_hess_for_irrep(
int irrep,
322 virtual void compute_mole(
const Displacement& d) =0;
327 unsigned int coor_to_irrep(
unsigned int symm_coord)
const;
333 class GradientsImpl :
public Impl {
345 int ndisplace()
const;
346 void set_gradient(
const Displacement& d,
double energy,
const RefSCVector &grad);
348 void compute_mole(
const Displacement& d);
354 class EnergiesImpl :
public Impl {
366 int ndisplace()
const;
368 void compute_mole(
const Displacement& d);
372 Eij(
int i,
int j, EnergiesImpl& eimpl) :
373 i_(i), j_(j), eimpl_(eimpl)
379 double operator()(
int di,
int dj) {
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)));
386 if (di+dj != 0) disp.push_back(std::make_pair(i_,
double(di+dj)));
388 eimpl_.compute_mole(disp);
389 const double energy = eimpl_.values().find(disp).second.
energy();
394 EnergiesImpl& eimpl_;
405 void override_default_params();
494 const Ref<Params>& params()
const {
return params_; }
529 double energy_accuracy_;
533 std::string restart_file_;
537 std::string checkpoint_file_;
540 int eliminate_quadratic_terms_;
546 std::vector<double> energies_;
549 void get_disp(
int disp,
int &index,
double &dispsize);
559 int ndisplace()
const;
560 int ndisplacements_done()
const {
return energies_.size(); }
562 void displace(
int disp);
563 void original_geometry();
565 void checkpoint_displacements(
StateOut&);
566 void restore_displacements(
StateIn&);
635 RefSCDimension d3natom()
const {
return mole_->moldim(); }
639 void set_eliminate_quadratic_terms(
bool e) { eliminate_quadratic_terms_ = e; }
640 void set_disp_size(
double s);