Psi4
optking/molecule.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * Psi4: an open-source quantum chemistry software package
5  *
6  * Copyright (c) 2007-2017 The Psi4 Developers.
7  *
8  * The copyrights for code used from other parties are included in
9  * the corresponding files.
10  *
11  * This program is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation; either version 2 of the License, or
14  * (at your option) any later version.
15  *
16  * This program 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 General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License along
22  * with this program; if not, write to the Free Software Foundation, Inc.,
23  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24  *
25  * @END LICENSE
26  */
27 
32 #ifndef _opt_molecule_h_
33 #define _opt_molecule_h_
34 
35 #include "package.h"
36 
37 #include "frag.h"
38 #include "interfrag.h"
39 #include "fb_frag.h"
40 #include "print.h"
42 #include <fstream>
43 
44 namespace opt {
45 
46 class MOLECULE {
47 
48  vector<FRAG *> fragments; // fragments with intrafragment coordinates
49  vector<INTERFRAG *> interfragments; // interfragment coordinates
50  vector<FB_FRAG *> fb_fragments; // EFP fixed-body fragment - does not actually
51  // store atoms (for now at least)
52  double energy;
53 
54  public:
55 
56  MOLECULE(int num_atoms); // allocate molecule with one fragment of this size
57 
59  for (std::size_t i=0; i<fragments.size(); ++i)
60  delete fragments[i];
61  fragments.clear();
62  for (std::size_t i=0; i<interfragments.size(); ++i)
63  delete interfragments[i];
64  interfragments.clear();
65  for (std::size_t i=0; i<fb_fragments.size(); ++i)
66  delete fb_fragments[i];
67  fb_fragments.clear();
68  }
69 
70  // if you have one fragment - make sure all atoms are bonded
71  // if not, break up into multiple fragments
72  void fragmentize(void);
73 
74  void add_interfragment(void);
75 
76  // Determine trivial, simple coordinate combinations.
78 
79  // Determine initial delocalized coordinate coefficients.
81 
82  // Determine Pulay natural coordinate combinations.
84 
85  int add_cartesians(void);
86 
87  int g_nfragment(void) const { return fragments.size(); };
88 
89  int g_nfb_fragment(void) const { return fb_fragments.size(); };
90 
91  int g_natom(void) const { // excludes atoms in fb fragments
92  int n = 0;
93  for (std::size_t f=0; f<fragments.size(); ++f)
94  n += fragments[f]->g_natom();
95  return n;
96  }
97 
98  int Ncoord(void) const {
99  int n=0;
100  for (std::size_t f=0; f<fragments.size(); ++f)
101  n += fragments[f]->Ncoord();
102  for (std::size_t i=0; i<interfragments.size(); ++i)
103  n += interfragments[i]->Ncoord();
104  for (std::size_t e=0; e<fb_fragments.size(); ++e)
105  n += fb_fragments[e]->Ncoord();
106  return n;
107  }
108 
109  int Ncoord_intrafragment(void) const {
110  int n=0;
111  for (std::size_t f=0; f<fragments.size(); ++f)
112  n += fragments[f]->Ncoord();
113  return n;
114  }
115 
116  int Ncoord_interfragment(void) const {
117  int n=0;
118  for (std::size_t f=0; f<interfragments.size(); ++f)
119  n += interfragments[f]->Ncoord();
120  return n;
121  }
122 
123  int Ncoord_fb_fragment(void) const {
124  int n=0;
125  for (std::size_t f=0; f<fb_fragments.size(); ++f)
126  n += fb_fragments[f]->Ncoord();
127  return n;
128  }
129 
130  // given fragment index returns first atom in that fragment
131  int g_atom_offset(int index) const {
132  int n = 0;
133  for (int f=1; f<=index; ++f)
134  n += fragments[f-1]->g_natom();
135  return n;
136  }
137 
138  // given fragment index, returns the absolute index
139  // for the first internal coordinate on that fragment
140  int g_coord_offset(int index) const {
141  int n = 0;
142  for (int f=0; f<index; ++f)
143  n += fragments[f]->Ncoord();
144  return n;
145  }
146 
147  // given interfragment index, returns the absolute index for the first
148  // internal coordinate of that interfragment set
149  int g_interfragment_coord_offset(int index) const {
150  int n = Ncoord_intrafragment();
151  for (int f=0; f<index; ++f)
152  n += interfragments[f]->Ncoord();
153  return n;
154  }
155 
156  // given fb_fragment index, returns the absolute index for the first
157  // internal coordinate of that fb set
158  int g_fb_fragment_coord_offset(int index) const {
160  for (int f=0; f<index; ++f)
161  n += fb_fragments[f]->Ncoord();
162  return n;
163  }
164 
165  double g_energy(void) const { return energy; }
166 
168  for (std::size_t i=0; i<fragments.size(); ++i)
170  }
171 
173  for (std::size_t i=0; i<fragments.size(); ++i)
175  }
176 
177  void print_connectivity(std::string psi_fp, FILE *qc_fp) const {
178  for (std::size_t i=0; i<fragments.size(); ++i)
179  fragments[i]->print_connectivity(psi_fp, qc_fp, i, g_atom_offset(i));
180  }
181 
182  void print_geom(std::string psi_fp, FILE *qc_fp, bool print_mass = false) {
183  for (std::size_t i=0; i<fragments.size(); ++i)
184  fragments[i]->print_geom(psi_fp, qc_fp, i, print_mass);
185  }
186 
187  void print_xyz(int iter_shift = 0);
188  void print_xyz_irc(int point, bool direction);
189 
190  void print_geom_grad(std::string psi_fp, FILE *qc_fp, bool print_mass = false) {
191  for (std::size_t i=0; i<fragments.size(); ++i)
192  fragments[i]->print_geom_grad(psi_fp, qc_fp, i, print_mass);
193  }
194 
195  void print_coords(std::string psi_fp, FILE *qc_fp) const;
196  void print_simples(std::string psi_fp, FILE *qc_fp) const;
197 
198  // print definition of an internal coordinate from global index
199  std::string get_coord_definition_from_global_index(int coord_index) const;
200 
201  void update_fb_values(void);
202 
203  void print_intco_dat(std::string psi_fp_coord, FILE *qc_fp);
204 
206  int n=0;
207  for (std::size_t i=0; i<fragments.size(); ++i)
208  n += fragments[i]->add_simples_by_connectivity();
209  return n;
210  }
211 
213  int n=0;
214  for (std::size_t i=0; i<fragments.size(); ++i)
215  n += fragments[i]->add_hbonds();
216  return n;
217  }
218 
220  int n=0;
221  for (std::size_t i=0; i<fragments.size(); ++i)
222  n += fragments[i]->add_auxiliary_bonds();
223  return n;
224  }
225 
226  // compute coord values from frag member geometries
227  double * coord_values(void) const {
228  GeomType x = g_geom_2D();
229  double *q = coord_values(x);
230  return q;
231  }
232 
233  // compute coord values from given geometry ; empty space for EFP coordinates included
234  double * coord_values(GeomType new_geom) const {
235  double *q, *q_frag, *q_IF;
236  q = init_array(Ncoord());
237 
238  for (std::size_t f=0; f<fragments.size(); ++f) {
239  q_frag = fragments[f]->coord_values( &(new_geom[g_atom_offset(f)]) );
240 
241  for (int i=0; i<fragments[f]->Ncoord(); ++i)
242  q[g_coord_offset(f)+i] = q_frag[i];
243 
244  free_array(q_frag);
245  }
246 
247  for (std::size_t I=0; I<interfragments.size(); ++I) {
248  int A_index = interfragments[I]->g_A_index();
249  int B_index = interfragments[I]->g_B_index();
250 
251  q_IF = interfragments[I]->coord_values( &(new_geom[g_atom_offset(A_index)]),
252  &(new_geom[g_atom_offset(B_index)]) );
253 
254  for (int i=0; i<interfragments[I]->Ncoord(); ++i)
255  q[g_interfragment_coord_offset(I)+i] = q_IF[i];
256 
257  free_array(q_IF);
258  }
259 
260  return q;
261  }
262 
263  void write_geom(void);
264  void symmetrize_geom(bool flexible=false);
265  void print_geom_out(void);
266  void print_geom_out_irc(void);
267 
268  double ** compute_B(void) const;
269  double ** compute_derivative_B(int coord_index) const ;
270 
271  double ** compute_G(bool use_masses=false) const;
272 
273  double * g_grad_array(void) const {
274  double *g, *g_frag;
275 
276  g = init_array(3*g_natom());
277  for (std::size_t f=0; f<fragments.size(); ++f) {
278  g_frag = fragments[f]->g_grad_array();
279  for (int i=0; i<3*fragments[f]->g_natom(); ++i)
280  g[3*g_atom_offset(f)+i] = g_frag[i];
281  free_array(g_frag);
282  }
283  return g;
284  }
285 
286  double * g_masses(void) const;
287  double * g_Z(void) const;
288  double * g_u_vector(void) const; // reciprocal masses in vector
289 
290  double * g_geom_array(void) {
291  double *g, *g_frag;
292 
293  g = init_array(3*g_natom());
294  for (std::size_t f=0; f<fragments.size(); ++f) {
295  g_frag = fragments[f]->g_geom_array();
296  for (int i=0; i<3*fragments[f]->g_natom(); ++i)
297  g[3*g_atom_offset(f)+i] = g_frag[i];
298  free_array(g_frag);
299  }
300  return g;
301  }
302 
303  double ** g_geom_2D(void) const {
304  double **g_frag;
305  double **g = init_matrix(g_natom(),3);
306 
307  for (std::size_t f=0; f<fragments.size(); ++f) {
308  g_frag = fragments[f]->g_geom();
309  for (int i=0; i<fragments[f]->g_natom(); ++i)
310  for (int xyz=0; xyz<3; ++xyz)
311  g[g_atom_offset(f)+i][xyz] = g_frag[i][xyz];
312  free_matrix(g_frag);
313  }
314  return g;
315  }
316 
317  double ** g_grad_2D(void) const {
318  double **g, *g_frag;
319 
320  g = init_matrix(g_natom(),3);
321  for (std::size_t f=0; f<fragments.size(); ++f) {
322  g_frag = fragments[f]->g_grad_array();
323  int cnt=0;
324  for (int i=0; i<fragments[f]->g_natom(); ++i)
325  for (int xyz=0; xyz<3; ++xyz)
326  g[g_atom_offset(f)+i][xyz] = g_frag[cnt++];
327  free_array(g_frag);
328  }
329  return g;
330  }
331 
332  void H_guess(void) const;
333  double **Lindh_guess(void) const;
334 
335  void forces(void);
336  void apply_constraint_forces(void);
337  bool has_fixed_eq_vals(void);
338  void project_f_and_H(void);
339  void project_dq(double *);
340  void irc_step(void);
341  void nr_step(void);
342  void rfo_step(void);
343  void prfo_step(void);
344  void backstep(void);
345  void sd_step(void);
346  //void sd_step_cartesians(void); now obsolete
347  void linesearch_step(void);
348 
349  void apply_intrafragment_step_limit(double * & dq);
350  std::vector<int> validate_angles(double const * const dq);
351 
352  void set_geom_array(double * array_in) {
353  for (std::size_t f=0; f<fragments.size(); ++f)
354  fragments[f]->set_geom_array( &(array_in[3*g_atom_offset(f)]) );
355  }
356 
357  void fix_tors_near_180(void) {
358  for (std::size_t f=0; f<fragments.size(); ++f)
360  for (std::size_t I=0; I<interfragments.size(); ++I)
362  }
363 
364  void fix_bend_axes(void) {
365  for (std::size_t f=0; f<fragments.size(); ++f)
366  fragments[f]->fix_bend_axes();
367  for (std::size_t I=0; I<interfragments.size(); ++I)
369  }
370  void unfix_bend_axes(void) {
371  for (std::size_t f=0; f<fragments.size(); ++f)
373  for (std::size_t I=0; I<interfragments.size(); ++I)
375  }
376 
377  void fix_oofp_near_180(void) {
378  for (std::size_t f=0; f<fragments.size(); ++f)
380  for (std::size_t I=0; I<interfragments.size(); ++I)
382  }
383 
384 /*
385  void check_tors_for_bad_angles(void) {
386  for (std::size_t f=0; f<fragments.size(); ++f)
387  fragments[f]->check_tors_for_bad_angles();
388  for (std::size_t I=0; I<interfragments.size(); ++I)
389  interfragments[I]->check_tors_for_bad_angles();
390  }
391 */
392 
393  void test_B(void);
394  void test_derivative_B(void);
395 
396  bool cartesian_H_to_internals(double **H_cart) const;
397 
398  void set_masses(void) {
399  for (std::size_t f=0; f<fragments.size(); ++f)
400  fragments[f]->set_masses();
401  }
402 
403  bool read_coords(std::ifstream & fin);
404  // function to obtain geometry and gradient
405  void read_geom_grad(void);
406 
407  // Compute constraint matrix
408  double ** compute_constraints(void);
409 
410  void add_fb_fragments(void);
411 
412  // freeze interfragment modes that break symmetry
413  void freeze_interfragment_asymm(void);
414 
415  // freeze all fragments in molecule
416  void freeze_intrafragments(void);
417 
418  // Freeze all coordinates within fragments
419  void freeze_intrafragment_coords(void);
420 
421  // determine whether a linear combination of coords breaks symmetry
422  bool coord_combo_is_symmetric(double *coord_combo, int dim);
423 
424  // Apply string list of user-specified internals to be frozen.
426 
427  // Tell if coord i is a fixed coordinate
428  bool is_coord_fixed(int coord_index);
429 
430  // Does the current set of coordinates includes ANYTHING that is not a cartesian.
431  bool is_noncart_present() const;
432 
433 };
434 
435 }
436 
437 #endif
double ** compute_G(bool use_masses=false) const
Definition: optking/molecule.cc:711
double ** compute_constraints(void)
Definition: optking/molecule.cc:866
double ** compute_derivative_B(int coord_index) const
Definition: optking/molecule.cc:607
const double *const *const GeomType
Definition: simple_base.h:44
double ** g_geom_2D(void) const
Definition: optking/molecule.h:303
int Ncoord_intrafragment(void) const
Definition: optking/molecule.h:109
void irc_step(void)
Definition: molecule_irc_step.cc:180
double * g_geom_array(void)
Definition: optking/molecule.h:290
int Ncoord_interfragment(void) const
Definition: optking/molecule.h:116
int g_natom(void) const
Definition: optking/molecule.h:91
double * g_masses(void) const
Definition: optking/molecule.cc:570
int add_intrafragment_hbonds(void)
Definition: optking/molecule.h:212
void update_fb_values(void)
Definition: optking/molecule.cc:985
bool read_coords(std::ifstream &fin)
Definition: molecule_read_coords.cc:116
void fix_tors_near_180(void)
Definition: optking/molecule.h:357
int Ncoord_fb_fragment(void) const
Definition: optking/molecule.h:123
void unfix_bend_axes(void)
Definition: optking/molecule.h:370
void backstep(void)
Definition: molecule_backstep.cc:89
Definition: pointgrp.h:105
bool coord_combo_is_symmetric(double *coord_combo, int dim)
Definition: optking/molecule.cc:902
void update_connectivity_by_distances(void)
Definition: optking/molecule.h:167
void read_geom_grad(void)
Definition: geom_gradients_io.cc:97
void fragmentize(void)
Definition: molecule_fragments.cc:68
int g_nfb_fragment(void) const
Definition: optking/molecule.h:89
void add_interfragment(void)
Definition: molecule_fragments.cc:290
vector< INTERFRAG * > interfragments
Definition: optking/molecule.h:49
int form_delocalized_coord_combinations(void)
Definition: optking/molecule.cc:796
int g_interfragment_coord_offset(int index) const
Definition: optking/molecule.h:149
Definition: optking/molecule.h:46
void linesearch_step(void)
Definition: molecule_linesearch_step.cc:64
void add_fb_fragments(void)
Definition: optking/molecule.cc:948
double * coord_values(void) const
Definition: optking/molecule.h:227
bool has_fixed_eq_vals(void)
Definition: optking/molecule.cc:157
void print_geom_out_irc(void)
Definition: molecule_print.cc:69
double ** compute_B(void) const
Definition: optking/molecule.cc:591
void print_xyz_irc(int point, bool direction)
Definition: molecule_print.cc:98
int g_nfragment(void) const
Definition: optking/molecule.h:87
void print_geom(std::string psi_fp, FILE *qc_fp, bool print_mass=false)
Definition: optking/molecule.h:182
int g_coord_offset(int index) const
Definition: optking/molecule.h:140
bool is_noncart_present() const
Definition: optking/molecule.cc:167
void nr_step(void)
Definition: molecule_nr_step.cc:59
void sd_step(void)
Definition: molecule_sd_step.cc:59
void write_geom(void)
Definition: geom_gradients_io.cc:341
void free_matrix(double **A)
Definition: mem.cc:137
void rfo_step(void)
Definition: molecule_rfo_step.cc:60
int form_natural_coord_combinations(void)
Definition: optking/molecule.cc:855
int form_trivial_coord_combinations(void)
Definition: optking/molecule.cc:786
void symmetrize_geom(bool flexible=false)
Definition: geom_gradients_io.cc:297
double * g_Z(void) const
Definition: optking/molecule.cc:579
void forces(void)
Definition: optking/molecule.cc:89
void project_f_and_H(void)
Definition: optking/molecule.cc:218
void print_connectivity(std::string psi_fp, FILE *qc_fp) const
Definition: optking/molecule.h:177
double ** init_matrix(long int m, long int n)
Definition: mem.cc:102
void fix_bend_axes(void)
Definition: optking/molecule.h:364
void freeze_intrafragment_coords(void)
Definition: optking/molecule.cc:779
int g_atom_offset(int index) const
Definition: optking/molecule.h:131
int add_intrafragment_simples_by_connectivity(void)
Definition: optking/molecule.h:205
void test_derivative_B(void)
Definition: molecule_tests.cc:151
int add_intrafragment_auxiliary_bonds(void)
Definition: optking/molecule.h:219
void free_array(double *f)
Definition: mem.cc:88
void apply_constraint_forces(void)
Definition: optking/molecule.cc:180
bool is_coord_fixed(int coord_index)
Definition: optking/molecule.cc:935
void print_intco_dat(std::string psi_fp_coord, FILE *qc_fp)
Definition: molecule_print.cc:129
~MOLECULE()
Definition: optking/molecule.h:58
void project_dq(double *)
Definition: optking/molecule.cc:316
int g_fb_fragment_coord_offset(int index) const
Definition: optking/molecule.h:158
void print_xyz(int iter_shift=0)
Definition: molecule_print.cc:82
void set_masses(void)
Definition: optking/molecule.h:398
void prfo_step(void)
Definition: molecule_prfo_step.cc:55
void set_geom_array(double *array_in)
Definition: optking/molecule.h:352
void freeze_intrafragments(void)
Definition: optking/molecule.cc:772
vector< FB_FRAG * > fb_fragments
Definition: optking/molecule.h:50
double energy
Definition: optking/molecule.h:52
double * g_u_vector(void) const
Definition: optking/molecule.cc:73
double * init_array(long int m)
Definition: mem.cc:59
void H_guess(void) const
Definition: optking/molecule.cc:415
MOLECULE(int num_atoms)
Definition: optking/molecule.cc:62
bool apply_input_constraints()
Definition: optking/molecule.cc:739
void print_simples(std::string psi_fp, FILE *qc_fp) const
Definition: molecule_print.cc:232
int add_cartesians(void)
Definition: optking/molecule.cc:764
std::string get_coord_definition_from_global_index(int coord_index) const
Definition: molecule_print.cc:153
void freeze_interfragment_asymm(void)
Definition: molecule_fragments.cc:532
void test_B(void)
Definition: molecule_tests.cc:58
std::vector< int > validate_angles(double const *const dq)
Definition: optking/molecule.cc:393
vector< FRAG * > fragments
Definition: optking/molecule.h:48
void apply_intrafragment_step_limit(double *&dq)
Definition: optking/molecule.cc:371
double * g_grad_array(void) const
Definition: optking/molecule.h:273
void fix_oofp_near_180(void)
Definition: optking/molecule.h:377
void print_coords(std::string psi_fp, FILE *qc_fp) const
Definition: molecule_print.cc:205
double g_energy(void) const
Definition: optking/molecule.h:165
double ** Lindh_guess(void) const
Definition: lindh_guess.cc:77
bool cartesian_H_to_internals(double **H_cart) const
Definition: optking/molecule.cc:481
double ** g_grad_2D(void) const
Definition: optking/molecule.h:317
void update_connectivity_by_bonds(void)
Definition: optking/molecule.h:172
double * coord_values(GeomType new_geom) const
Definition: optking/molecule.h:234
void print_geom_grad(std::string psi_fp, FILE *qc_fp, bool print_mass=false)
Definition: optking/molecule.h:190
int Ncoord(void) const
Definition: optking/molecule.h:98
void print_geom_out(void)
Definition: molecule_print.cc:58