1 // 2 // fdhess.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_fdhess_h 29 #define _chemistry_molecule_fdhess_h 30 31 #ifdef __GNUC__ 32 #pragma interface 33 #endif 34 35 #include <iostream> 36 37 #include <chemistry/molecule/hess.h> 38 #include <chemistry/molecule/energy.h> 39 40 namespace sc { 41 42 /** Computes the molecular hessian by finite displacements of gradients. 43 This will use the minimum number of displacements, each in the 44 highest possible point group. */ 45 class FinDispMolecularHessian: public MolecularHessian { 46 protected: 47 Ref<MolecularEnergy> mole_; 48 // In case molecule must be given in lower symmetry, its actual 49 // symmetry and the symmetry used to compute displacements is this 50 Ref<PointGroup> displacement_point_group_; 51 // The molecule's original point group for restoration at the end. 52 Ref<PointGroup> original_point_group_; 53 // The molecule's original geometry for restoration at the end and 54 //computing displacements. 55 RefSCVector original_geometry_; 56 // the cartesian displacement size in bohr 57 double disp_; 58 // the accuracy for gradient calculations 59 double accuracy_; 60 // the number of completed displacements 61 int ndisp_; 62 // the number of irreps in the displacement point group 63 int nirrep_; 64 // whether or not to attempt a restart 65 int restart_; 66 // the name of the restart file 67 char *restart_file_; 68 // whether or not to checkpoint 69 int checkpoint_; 70 // the name of the checkpoint file 71 char *checkpoint_file_; 72 // only do the totally symmetric displacements 73 int only_totally_symmetric_; 74 // eliminate the cubic terms by doing an extra displacement for 75 //each of the totally symmetry coordinates 76 int eliminate_cubic_terms_; 77 // use the gradient at the initial geometry to remove first order terms 78 // (important if not at equilibrium geometry) 79 int do_null_displacement_; 80 // print flag 81 int debug_; 82 // a basis for the symmetrized cartesian coordinates 83 RefSCMatrix symbasis_; 84 // the gradients at each of the displacements 85 RefSCVector *gradients_; 86 87 void get_disp(int disp, int &irrep, int &index, double &coef); 88 void do_hess_for_irrep(int irrep, 89 const RefSymmSCMatrix &dhessian, 90 const RefSymmSCMatrix &xhessian); 91 void init(); 92 void restart(); 93 public: 94 FinDispMolecularHessian(const Ref<MolecularEnergy>&); 95 /** The FinDispMolecularHessian KeyVal constructor is used to generate a 96 FinDispMolecularHessian object from the input. It reads the keywords 97 below. 98 99 <table border="1"> 100 101 <tr><td>Keyword<td>Type<td>Default<td>Description 102 103 <tr><td><tt>energy</tt><td>MolecularEnergy<td>none<td>This gives an 104 object which will be used to compute the gradients needed to form 105 the hessian. If this is not specified, the object using 106 FinDispMolecularHessian will, in some cases, fill it in 107 appropriately. However, even in these cases, it may be desirable 108 to specify this keyword. For example, this could be used in an 109 optimization to compute frequencies using a lower level of theory. 110 111 <tr><td><tt>debug</tt><td>boolean<td>false<td>If true, 112 print out debugging information. 113 114 <tr><td><tt>point_group</tt><td>PointGroup<td>none<td> 115 The point group to use for generating the displacements. 116 117 <tr><td><tt>restart</tt><td>boolean<td>true<td>If true, and a 118 checkpoint file exists, restart from that file. 119 120 <tr><td><tt>restart_file</tt><td>string 121 <td><em>basename</em><tt>.ckpt.hess</tt><td>The name of 122 the file where checkpoint information is written to or read from. 123 124 <tr><td><tt>checkpoint</tt><td>boolean<td>true<td>If true, 125 checkpoint intermediate data. 126 127 <tr><td><tt>only_totally_symmetric</tt><td>boolean<td>false 128 <td>If true, only follow totally symmetric displacments. The 129 hessian will not be complete, but it has enough information 130 to use it in a geometry optimization. 131 132 <tr><td><tt>eliminate_cubic_terms</tt><td>boolean<td>true<td> 133 If true, then cubic terms will be eliminated. This requires 134 that two displacements are done for each totally symmetric 135 coordinate, rather than one. Setting this to false will reduce 136 the accuracy, but the results will still probably be accurate 137 enough for a geometry optimization. 138 139 <tr><td><tt>do_null_displacement</tt><td>boolean<td>true<td>Run 140 the calculation at the given geometry as well. 141 142 <tr><td><tt>displacement</tt><td>double<td>1.0e-2<td>The size of 143 the displacement in Bohr. 144 145 <tr><td><tt>gradient_accuracy</tt><td>double<td><tt>displacement</tt> 146 / 1000<td>The accuracy to which the gradients will be computed. 147 148 </table> 149 */ 150 FinDispMolecularHessian(const Ref<KeyVal>&); 151 FinDispMolecularHessian(StateIn&); 152 ~FinDispMolecularHessian(); 153 void save_data_state(StateOut&); 154 155 /** These members are used to compute a cartesian hessian from 156 gradients at finite displacements. */ 157 RefSymmSCMatrix compute_hessian_from_gradients(); 158 int ndisplace() const; ndisplacements_done()159 int ndisplacements_done() const { return ndisp_; } 160 RefSCMatrix displacements(int irrep) const; 161 void displace(int disp); 162 void original_geometry(); 163 void set_gradient(int disp, const RefSCVector &grad); 164 void checkpoint_displacements(StateOut&); 165 void restore_displacements(StateIn&); 166 167 /** This returns the cartesian hessian. If it has not yet been 168 computed, it will be computed by finite displacements. */ 169 RefSymmSCMatrix cartesian_hessian(); 170 171 /// Set checkpoint option. set_checkpoint(int c)172 void set_checkpoint(int c) { checkpoint_ = c; } 173 /// Return the current value of the checkpoint option. checkpoint()174 int checkpoint() const { return checkpoint_; } 175 176 void set_energy(const Ref<MolecularEnergy> &energy); 177 MolecularEnergy* energy() const; 178 matrixkit()179 Ref<SCMatrixKit> matrixkit() const { return mole_->matrixkit(); } d3natom()180 RefSCDimension d3natom() const { return mole_->moldim(); } 181 }; 182 183 } 184 185 #endif 186 187 // Local Variables: 188 // mode: c++ 189 // c-file-style: "CLJ" 190 // End: 191