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:
int equivalent(int iuniq, int j) const
Returns the j'th atom equivalent to iuniq.
Definition: molecule.h:406
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:433
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:404
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:402
bool include_qq() const
Returns include_qq. See set_include_qq.
Definition: molecule.h:439
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:431
int q_atom(int i) const
Retrieve the "Q" atoms.
Definition: molecule.h:444
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:409
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:449
int nunique() const
Return information about symmetry unique and equivalent atoms.
Definition: molecule.h:400
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:437
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.
SCVector3 center_of_charge() const
Returns a SCVector3 containing the cartesian coordinates of the center of charge for 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:447
Ref< AtomInfo > atominfo() const
Return the molecule's AtomInfo object.
Definition: molecule.h:421
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:442

Generated at Mon Jun 27 2016 11:47:20 for MPQC 2.3.1 using the documentation package Doxygen 1.8.10.