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 
354  double nuclear_repulsion_energy();
355 
358  void nuclear_repulsion_1der(int center, double xyz[3]);
359 
361  void nuclear_efield(const double *position, double* efield);
362 
365  void nuclear_charge_efield(const double *charges,
366  const double *position, double* efield);
367 
373  void symmetrize(double tol = 0.5);
374 
376  void symmetrize(const Ref<PointGroup> &pg, double tol = 0.5);
377 
381  void cleanup_molecule(double tol = 0.1);
382 
383  void translate(const double *r);
384  void move_to_com();
385  void transform_to_principal_axes(int trans_frame=1);
386  void transform_to_symmetry_frame();
387  void print_pdb(std::ostream& =ExEnv::out0(), char *title =0) const;
388 
389  void read_pdb(const char *filename);
390 
393  void principal_moments_of_inertia(double *evals, double **evecs=0) const;
394 
396  int nunique() const { return nuniq_; }
398  int unique(int iuniq) const { return equiv_[iuniq][0]; }
400  int nequivalent(int iuniq) const { return nequiv_[iuniq]; }
402  int equivalent(int iuniq, int j) const { return equiv_[iuniq][j]; }
405  int atom_to_unique(int iatom) const { return atom_to_uniq_[iatom]; }
408  int atom_to_unique_offset(int iatom) const;
409 
411  int n_core_electrons();
412 
414  int max_z();
415 
417  Ref<AtomInfo> atominfo() const { return atominfo_; }
418 
420  std::string atom_name(int iatom) const;
421 
423  std::string atom_symbol(int iatom) const;
424 
427  void set_include_q(bool iq) { include_q_ = iq; }
429  bool include_q() const { return include_q_; }
430 
433  void set_include_qq(bool iqq) { include_qq_ = iqq; }
435  bool include_qq() const { return include_qq_; }
436 
438  int n_q_atom() const { return q_atoms_.size(); }
440  int q_atom(int i) const { return q_atoms_[i]; }
441 
443  int n_non_q_atom() const { return non_q_atoms_.size(); }
445  int non_q_atom(int i) const { return non_q_atoms_[i]; }
446 
447  void save_data_state(StateOut&);
448 };
449 
450 }
451 
452 #endif
453 
454 // Local Variables:
455 // mode: c++
456 // c-file-style: "CLJ"
457 // End:
int equivalent(int iuniq, int j) const
Returns the j'th atom equivalent to iuniq.
Definition: molecule.h:402
void nuclear_repulsion_1der(int center, double xyz[3])
Compute the nuclear repulsion energy first derivative with respect to the given center.
void nuclear_efield(const double *position, double *efield)
Compute the electric field due to the nuclei at the given point.
Ref< PointGroup > highest_point_group(double tol=1.0e-8) const
Find this molecules true point group (limited to abelian groups).
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
std::string atom_name(int iatom) const
Returns the element name of the atom.
int is_plane(SCVector3 &origin, SCVector3 &uperp, double tol=1.0e-8) const
Return 1 if the given plane is a symmetry element for the molecule.
int is_planar(double tolerance=1.0e-5) const
Returns 1 if the molecule is planar, 0 otherwise.
void is_linear_planar(int &linear, int &planar, double tol=1.0e-5) const
Sets linear to 1 if the molecular is linear, 0 otherwise.
bool include_q() const
Returns include_q. See set_include_q.
Definition: molecule.h:429
Serializes objects that derive from SavableState.
Definition: stateout.h:61
int nequivalent(int iuniq) const
Returns the number of atoms equivalent to iuniq.
Definition: molecule.h:400
void nuclear_charge_efield(const double *charges, const double *position, double *efield)
Compute the electric field due to the given charges at the positions of the nuclei at the given point...
int unique(int iuniq) const
Returns the overall number of the iuniq'th unique atom.
Definition: molecule.h:398
bool include_qq() const
Returns include_qq. See set_include_qq.
Definition: molecule.h:435
void set_include_q(bool iq)
If include_q is true, then include the "Q" atoms in the charge and efield routines.
Definition: molecule.h:427
int q_atom(int i) const
Retrieve the "Q" atoms.
Definition: molecule.h:440
A template class that maintains references counts.
Definition: ref.h:332
double nuclear_repulsion_energy()
Returns the nuclear repulsion energy for the molecule.
int is_linear(double tolerance=1.0e-5) const
Returns 1 if the molecule is linear, 0 otherwise.
int atom_to_unique(int iatom) const
Converts an atom number to the number of its generating unique atom.
Definition: molecule.h:405
Ref< PointGroup > point_group() const
Returns the PointGroup of the molecule.
void principal_moments_of_inertia(double *evals, double **evecs=0) const
Compute the principal moments of inertia and, possibly, the principal axes.
int atom_to_unique_offset(int iatom) const
Converts an atom number to the offset of this atom in the list of generated atoms.
Definition: mpqcin.h:13
void cleanup_molecule(double tol=0.1)
This will try to carefully correct symmetry errors in molecules.
Restores objects that derive from SavableState.
Definition: statein.h:70
int has_inversion(SCVector3 &origin, double tol=1.0e-8) const
Return 1 if the molecule has an inversion center.
std::string atom_symbol(int iatom) const
Returns the element symbol of the atom.
void symmetrize(double tol=0.5)
If the molecule contains only symmetry unique atoms, this function will generate the other...
SCVector3 center_of_mass() const
Returns a SCVector3 containing the cartesian coordinates of the center of mass for the molecule...
int atom_label_to_index(const char *label) const
Returns the index of the atom with the given label.
static std::ostream & out0()
Return an ostream that writes from node 0.
int atom_at_position(double *, double tol=0.05) const
Takes an (x, y, z) postion and finds an atom within the given tolerance distance. ...
double nuclear_charge() const
Returns the total nuclear charge.
int non_q_atom(int i) const
Retrieve the of non-"Q" atoms.
Definition: molecule.h:445
int nunique() const
Return information about symmetry unique and equivalent atoms.
Definition: molecule.h:396
void set_include_qq(bool iqq)
If include_qq is true, include the coupling between pairs of "Q" atoms when computing nuclear repulsi...
Definition: molecule.h:433
void add_atom(int Z, double x, double y, double z, const char *=0, double mass=0.0, int have_charge=0, double charge=0.0)
Add an AtomicCenter to the Molecule.
int n_core_electrons()
Return the number of core electrons.
int natom() const
Returns the number of atoms in the molcule.
Definition: molecule.h:288
The Molecule class contains information about molecules.
Definition: molecule.h:127
int is_axis(SCVector3 &origin, SCVector3 &udirection, int order, double tol=1.0e-8) const
Return 1 if this given axis is a symmetry element for the molecule.
int n_non_q_atom() const
Retrieve the number of non-"Q" atoms.
Definition: molecule.h:443
Ref< AtomInfo > atominfo() const
Return the molecule's AtomInfo object.
Definition: molecule.h:417
double charge(int iatom) const
Return the charge of the atom.
double * charges() const
Returns a double* containing the nuclear charges of the atoms.
Base class for objects that can save/restore state.
Definition: state.h:46
int max_z()
Return the maximum atomic number.
void set_point_group(const Ref< PointGroup > &, double tol=1.0e-7)
Sets the PointGroup of the molecule.
const char * label(int atom) const
Returns the label explicitly assigned to atom.
virtual void print(std::ostream &=ExEnv::out0()) const
Print information about the molecule.
int n_q_atom() const
Retrieve the number of "Q" atoms.
Definition: molecule.h:438

Generated at Thu Sep 11 2014 22:18:27 for MPQC 2.3.1 using the documentation package Doxygen 1.8.8.