MPQC  2.3.1
molecule.h
1 //
2 // molecule.h
3 //
4 // Copyright (C) 1996 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_molecule_h
29 #define _chemistry_molecule_molecule_h
30 
31 #ifdef __GNUC__
32 #pragma interface
33 #endif
34 
35 #include <stdio.h>
36 #include <iostream>
37 #include <util/class/class.h>
38 #include <util/state/state.h>
39 #include <util/keyval/keyval.h>
40 #include <util/misc/units.h>
41 #include <math/symmetry/pointgrp.h>
42 #include <math/scmat/vector3.h>
43 #include <math/scmat/matrix.h>
44 #include <chemistry/molecule/atominfo.h>
45 
46 namespace sc {
47 
127 class Molecule: public SavableState
128 {
129  protected:
130  int natoms_;
131  Ref<AtomInfo> atominfo_;
132  Ref<PointGroup> pg_;
133  Ref<Units> geometry_units_;
134  double **r_;
135  int *Z_;
136  double *charges_;
137 
138  // symmetry equiv info
139  int nuniq_;
140  int *nequiv_;
141  int **equiv_;
142  int *atom_to_uniq_;
143  void init_symmetry_info(double tol=0.5);
144  void clear_symmetry_info();
145 
146  // these are optional
147  double *mass_;
148  char **labels_;
149 
150  // The Z that represents a "Q" type atom.
151  int q_Z_;
152 
153  // If true, include the q terms in the charge and efield routines
154  bool include_q_;
155 
156  // If true, include the coupling between q-q pairs when
157  // computing nuclear repulsion energy and gradients.
158  bool include_qq_;
159 
160  // These vectors contain the atom indices of atoms that are not type
161  // "Q" and those that are.
162  std::vector<int> q_atoms_;
163  std::vector<int> non_q_atoms_;
164 
165  void clear();
166 
167  // Throw an exception if an atom is duplicated. The
168  // atoms in the range [begin, natom_) are checked.
169  void throw_if_atom_duplicated(int begin=0, double tol = 1e-3);
170  public:
171  Molecule();
172  Molecule(const Molecule&);
173  Molecule(StateIn&);
269  Molecule(const Ref<KeyVal>&input);
270 
271  virtual ~Molecule();
272 
273  Molecule& operator=(const Molecule&);
274 
276  void add_atom(int Z,double x,double y,double z,
277  const char * = 0, double mass = 0.0,
278  int have_charge = 0, double charge = 0.0);
279 
281  virtual void print(std::ostream& =ExEnv::out0()) const;
282  virtual void print_parsedkeyval(std::ostream& =ExEnv::out0(),
283  int print_pg = 1,
284  int print_unit = 1,
285  int number_atoms = 1) const;
286 
288  int natom() const { return natoms_; }
289 
290  int Z(int atom) const { return Z_[atom]; }
291  double &r(int atom, int xyz) { return r_[atom][xyz]; }
292  const double &r(int atom, int xyz) const { return r_[atom][xyz]; }
293  double *r(int atom) { return r_[atom]; }
294  const double *r(int atom) const { return r_[atom]; }
295  double mass(int atom) const;
298  const char *label(int atom) const;
299 
302  int atom_at_position(double *, double tol = 0.05) const;
303 
306  int atom_label_to_index(const char *label) const;
307 
311  double *charges() const;
312 
314  double charge(int iatom) const;
315 
317  double nuclear_charge() const;
318 
320  void set_point_group(const Ref<PointGroup>&, double tol=1.0e-7);
322  Ref<PointGroup> point_group() const;
323 
327  Ref<PointGroup> highest_point_group(double tol = 1.0e-8) const;
328 
331  int is_axis(SCVector3 &origin,
332  SCVector3 &udirection, int order, double tol=1.0e-8) const;
333 
336  int is_plane(SCVector3 &origin, SCVector3 &uperp, double tol=1.0e-8) const;
337 
339  int has_inversion(SCVector3 &origin, double tol = 1.0e-8) const;
340 
342  int is_linear(double tolerance = 1.0e-5) const;
344  int is_planar(double tolerance = 1.0e-5) const;
347  void is_linear_planar(int&linear,int&planar,double tol = 1.0e-5) const;
348 
351  SCVector3 center_of_mass() const;
352 
355  SCVector3 center_of_charge() const;
356 
358  double nuclear_repulsion_energy();
359 
362  void nuclear_repulsion_1der(int center, double xyz[3]);
363 
365  void nuclear_efield(const double *position, double* efield);
366 
369  void nuclear_charge_efield(const double *charges,
370  const double *position, double* efield);
371 
377  void symmetrize(double tol = 0.5);
378 
380  void symmetrize(const Ref<PointGroup> &pg, double tol = 0.5);
381 
385  void cleanup_molecule(double tol = 0.1);
386 
387  void translate(const double *r);
388  void move_to_com();
389  void transform_to_principal_axes(int trans_frame=1);
390  void transform_to_symmetry_frame();
391  void print_pdb(std::ostream& =ExEnv::out0(), char *title =0) const;
392 
393  void read_pdb(const char *filename);
394 
397  void principal_moments_of_inertia(double *evals, double **evecs=0) const;
398 
400  int nunique() const { return nuniq_; }
402  int unique(int iuniq) const { return equiv_[iuniq][0]; }
404  int nequivalent(int iuniq) const { return nequiv_[iuniq]; }
406  int equivalent(int iuniq, int j) const { return equiv_[iuniq][j]; }
409  int atom_to_unique(int iatom) const { return atom_to_uniq_[iatom]; }
412  int atom_to_unique_offset(int iatom) const;
413 
415  int n_core_electrons();
416 
418  int max_z();
419 
421  Ref<AtomInfo> atominfo() const { return atominfo_; }
422 
424  std::string atom_name(int iatom) const;
425 
427  std::string atom_symbol(int iatom) const;
428 
431  void set_include_q(bool iq) { include_q_ = iq; }
433  bool include_q() const { return include_q_; }
434 
437  void set_include_qq(bool iqq) { include_qq_ = iqq; }
439  bool include_qq() const { return include_qq_; }
440 
442  int n_q_atom() const { return q_atoms_.size(); }
444  int q_atom(int i) const { return q_atoms_[i]; }
445 
447  int n_non_q_atom() const { return non_q_atoms_.size(); }
449  int non_q_atom(int i) const { return non_q_atoms_[i]; }
450 
451  void save_data_state(StateOut&);
452 };
453 
454 }
455 
456 #endif
457 
458 // Local Variables:
459 // mode: c++
460 // c-file-style: "CLJ"
461 // End:

Generated at Tue Aug 20 2013 22:07:00 for MPQC 2.3.1 using the documentation package Doxygen 1.8.3.1.