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