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-2018 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();
74 
75  void add_interfragment();
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();
87 
88  int g_nfragment() const { return fragments.size(); };
89 
90  int g_nfb_fragment() const { return fb_fragments.size(); };
91 
92  int g_natom() 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() 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() 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() 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() 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() 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();
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() 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();
265  void symmetrize_geom(bool flexible=false);
266  void print_geom_out();
267  void print_geom_out_irc();
268 
269  double ** compute_B() 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() 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() const;
288  double * g_Z() const;
289  double * g_u_vector() const; // reciprocal masses in vector
290 
291  double * g_geom_array() {
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() 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() 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() const;
334  double **Lindh_guess() const;
335 
336  void forces();
338  bool has_fixed_eq_vals();
339  void project_f_and_H();
340  void project_dq(double *);
341  void irc_step();
342  void nr_step();
343  void rfo_step();
344  void prfo_step();
345  void backstep();
346  void sd_step();
347  //void sd_step_cartesians(void); now obsolete
348  void linesearch_step();
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 
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() {
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  }
372  for (std::size_t f=0; f<fragments.size(); ++f)
374  for (std::size_t I=0; I<interfragments.size(); ++I)
376  }
377 
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();
395  void test_derivative_B();
396 
397  bool cartesian_H_to_internals(double **H_cart) const;
398 
399  void set_masses() {
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();
407 
408  // Compute constraint matrix
409  double ** compute_constraints();
410 
411  void add_fb_fragments();
412 
413  // freeze interfragment modes that break symmetry
415 
416  // freeze all fragments in molecule
417  void freeze_intrafragments();
418 
419  // Freeze all coordinates within fragments
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:711
void irc_step()
Definition: molecule_irc_step.cc:181
void prfo_step()
Definition: molecule_prfo_step.cc:56
double ** compute_derivative_B(int coord_index) const
Definition: optking/molecule.cc:607
const double *const *const GeomType
Definition: simple_base.h:45
void add_interfragment()
Definition: molecule_fragments.cc:292
int form_delocalized_coord_combinations()
Definition: optking/molecule.cc:796
void set_masses()
Definition: optking/molecule.h:399
int Ncoord_interfragment() const
Definition: optking/molecule.h:117
double * coord_values() const
Definition: optking/molecule.h:228
void unfix_bend_axes()
Definition: optking/molecule.h:371
void print_geom_out()
Definition: molecule_print.cc:57
void test_derivative_B()
Definition: molecule_tests.cc:150
bool read_coords(std::ifstream &fin)
Definition: molecule_read_coords.cc:115
double * g_grad_array() const
Definition: optking/molecule.h:274
void sd_step()
Definition: molecule_sd_step.cc:60
int add_cartesians()
Definition: optking/molecule.cc:764
void nr_step()
Definition: molecule_nr_step.cc:60
double * g_Z() const
Definition: optking/molecule.cc:579
double ** Lindh_guess() const
Definition: lindh_guess.cc:78
Definition: pointgrp.h:104
bool coord_combo_is_symmetric(double *coord_combo, int dim)
Definition: optking/molecule.cc:902
int g_natom() const
Definition: optking/molecule.h:92
void project_f_and_H()
Definition: optking/molecule.cc:218
void fix_bend_axes()
Definition: optking/molecule.h:365
vector< INTERFRAG * > interfragments
Definition: optking/molecule.h:50
double ** g_grad_2D() const
Definition: optking/molecule.h:318
double * g_u_vector() const
Definition: optking/molecule.cc:73
int g_interfragment_coord_offset(int index) const
Definition: optking/molecule.h:150
int Ncoord_intrafragment() const
Definition: optking/molecule.h:110
Definition: optking/molecule.h:47
void fragmentize()
Definition: molecule_fragments.cc:70
void forces()
Definition: optking/molecule.cc:89
void rfo_step()
Definition: molecule_rfo_step.cc:61
int form_trivial_coord_combinations()
Definition: optking/molecule.cc:786
void print_xyz_irc(int point, bool direction)
Definition: molecule_print.cc:97
void print_geom(std::string psi_fp, FILE *qc_fp, bool print_mass=false)
Definition: optking/molecule.h:183
void H_guess() const
Definition: optking/molecule.cc:415
int g_coord_offset(int index) const
Definition: optking/molecule.h:141
bool is_noncart_present() const
Definition: optking/molecule.cc:167
void backstep()
Definition: molecule_backstep.cc:90
double ** compute_constraints()
Definition: optking/molecule.cc:866
double g_energy() const
Definition: optking/molecule.h:166
void free_matrix(double **A)
Definition: mem.cc:138
long int size_t
Definition: sortintegrals.cc:41
int g_nfragment() const
Definition: optking/molecule.h:88
int Ncoord_fb_fragment() const
Definition: optking/molecule.h:124
void freeze_intrafragments()
Definition: optking/molecule.cc:772
void linesearch_step()
Definition: molecule_linesearch_step.cc:65
void symmetrize_geom(bool flexible=false)
Definition: geom_gradients_io.cc:296
int add_intrafragment_simples_by_connectivity()
Definition: optking/molecule.h:206
void write_geom()
Definition: geom_gradients_io.cc:340
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 apply_constraint_forces()
Definition: optking/molecule.cc:180
int g_atom_offset(int index) const
Definition: optking/molecule.h:132
double * g_masses() const
Definition: optking/molecule.cc:570
void free_array(double *f)
Definition: mem.cc:89
double ** compute_B() const
Definition: optking/molecule.cc:591
int Ncoord() const
Definition: optking/molecule.h:99
bool is_coord_fixed(int coord_index)
Definition: optking/molecule.cc:935
double ** g_geom_2D() const
Definition: optking/molecule.h:304
int g_nfb_fragment() const
Definition: optking/molecule.h:90
void print_geom_out_irc()
Definition: molecule_print.cc:68
void update_connectivity_by_distances()
Definition: optking/molecule.h:168
void print_intco_dat(std::string psi_fp_coord, FILE *qc_fp)
Definition: molecule_print.cc:128
~MOLECULE()
Definition: optking/molecule.h:59
void project_dq(double *)
Definition: optking/molecule.cc:316
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:81
void set_geom_array(double *array_in)
Definition: optking/molecule.h:353
vector< FB_FRAG * > fb_fragments
Definition: optking/molecule.h:51
double energy
Definition: optking/molecule.h:53
void add_fb_fragments()
Definition: optking/molecule.cc:948
double * g_geom_array()
Definition: optking/molecule.h:291
double * init_array(long int m)
Definition: mem.cc:60
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:231
std::string get_coord_definition_from_global_index(int coord_index) const
Definition: molecule_print.cc:152
bool has_fixed_eq_vals()
Definition: optking/molecule.cc:157
void fix_tors_near_180()
Definition: optking/molecule.h:358
void update_fb_values()
Definition: optking/molecule.cc:985
std::vector< int > validate_angles(double const *const dq)
Definition: optking/molecule.cc:393
void read_geom_grad()
Definition: geom_gradients_io.cc:96
int add_intrafragment_hbonds()
Definition: optking/molecule.h:213
vector< FRAG * > fragments
Definition: optking/molecule.h:49
void fix_oofp_near_180()
Definition: optking/molecule.h:378
void apply_intrafragment_step_limit(double *&dq)
Definition: optking/molecule.cc:371
int add_intrafragment_auxiliary_bonds()
Definition: optking/molecule.h:220
void freeze_intrafragment_coords()
Definition: optking/molecule.cc:779
void print_coords(std::string psi_fp, FILE *qc_fp) const
Definition: molecule_print.cc:204
bool cartesian_H_to_internals(double **H_cart) const
Definition: optking/molecule.cc:481
void freeze_interfragment_asymm()
Definition: molecule_fragments.cc:534
void update_connectivity_by_bonds()
Definition: optking/molecule.h:173
int form_natural_coord_combinations()
Definition: optking/molecule.cc:855
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
void test_B()
Definition: molecule_tests.cc:57