molecule.h

//
// molecule.h
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <cljanss@limitpt.com>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//
 
#ifndef _chemistry_molecule_molecule_h
#define _chemistry_molecule_molecule_h
 
#undef HAVE_OPENBABEL2  // Won't link on my machine.  Remove this before committing to main branch
 
#include <stdio.h>
#include <iostream>
#include <mpqc_config.h>
#include <util/class/class.h>
#include <util/state/state.h>
#include <util/keyval/keyval.h>
#include <util/misc/units.h>
#include <math/symmetry/pointgrp.h>
#include <math/scmat/vector3.h>
#include <math/scmat/matrix.h>
#include <chemistry/molecule/atominfo.h>
#include <chemistry/molecule/atom.h>
 
#ifdef MPQC_NEW_FEATURES
#include <util/misc/xml.h>
#endif // MPQC_NEW_FEATURES
 
namespace sc {
 
 
class Molecule: public SavableState
#ifdef MPQC_NEW_FEATURES
, virtual public DescribedXMLWritable
#endif
{
  protected:
    std::vector<Atom> atoms_;
    Ref<AtomInfo> atominfo_;
    Ref<PointGroup> pg_;
    Ref<Units> geometry_units_;
    double ref_origin_[3];   //< position of the origin of the reference coordinate system
 
    // symmetry equiv info
    int nuniq_;
    int *nequiv_;
    int **equiv_;
    int *atom_to_uniq_;
    void init_symmetry_info(double tol=0.5);
    void clear_symmetry_info();
 
    // The Z that represents a "Q" type atom.
    int q_Z_;
 
    // If true, include the q terms in the charge and efield routines
    bool include_q_;
 
    // If true, include the coupling between q-q pairs when
    // computing nuclear repulsion energy and gradients.
    bool include_qq_;
 
    // These vectors contain the atom indices of atoms that are not type
    // "Q" and those that are.
    std::vector<int> q_atoms_;
    std::vector<int> non_q_atoms_;
 
    void clear();
 
    // Throw an exception if an atom is duplicated.  The
    // atoms in the range [begin, natom_) are checked.
    void throw_if_atom_duplicated(int begin=0, double tol = 1e-3);
  public:
    Molecule();
    Molecule(const Molecule&);
    Molecule(StateIn&);
    Molecule(const Ref<KeyVal>&input);
 
    virtual ~Molecule();
 
    Molecule& operator=(const Molecule&);
 
    void add_atom(int Z,double x,double y,double z,
                  const std::string & label = "", double mass = 0.0,
                  int have_charge = 0, double charge = 0.0,
                  int have_fragment = 0, int fragment = 0);
 
    const Atom& atom(size_t i) const { return atoms_[i]; }
    const std::vector<Atom>& atoms() const { return atoms_; }
 
    virtual void print(std::ostream& =ExEnv::out0()) const;
    virtual void print_parsedkeyval(std::ostream& =ExEnv::out0(),
                                    int print_pg = 1,
                                    int print_unit = 1,
                                    int number_atoms = 1) const;
 
    Ref<Units> geometry_units() const { return geometry_units_; }
 
    size_t natom() const { return atoms_.size(); }
 
    int Z(int atom) const { return atoms_[atom].Z(); }
    double &r(int atom, int xyz) { return atoms_[atom].r(xyz); }
    const double &r(int atom, int xyz) const { return atoms_[atom].r(xyz); }
    const double *r(int atom) const { return atoms_[atom].r(); }
    double mass(int atom) const;
    const char *label(int atom) const;
    int fragment(int atom) const;
 
    int atom_at_position(double *, double tol = 0.05) const;
 
    int atom_label_to_index(const std::string &label) const;
 
    std::vector<double> charges() const;
 
    double charge(int iatom) const;
 
    bool is_Q(int iatom) const;
 
    double total_charge() const;
 
    int total_Z() const;
 
    void set_point_group(const Ref<PointGroup>&, double tol=1.0e-7);
    const Ref<PointGroup>& point_group() const;
 
    Ref<PointGroup> highest_point_group(double tol = 1.0e-8) const;
 
    int is_axis(SCVector3 &origin,
                SCVector3 &udirection, int order, double tol=1.0e-8) const;
 
    int is_plane(SCVector3 &origin, SCVector3 &uperp, double tol=1.0e-8) const;
 
    int has_inversion(SCVector3 &origin, double tol = 1.0e-8) const;
 
    int is_linear(double tolerance = 1.0e-5) const;
    int is_planar(double tolerance = 1.0e-5) const;
    void is_linear_planar(int&linear,int&planar,double tol = 1.0e-5) const;
 
    SCVector3 center_of_mass() const;
 
    double nuclear_repulsion_energy();
 
    void nuclear_repulsion_1der(int center, double xyz[3]);
 
    void nuclear_efield(const double *position, double* efield);
 
    void nuclear_charge_efield(const double *charges,
                               const double *position, double* efield);
 
    void symmetrize(double tol = 0.5);
 
    void symmetrize(const Ref<PointGroup> &pg, double tol = 0.5);
 
    void cleanup_molecule(double tol = 0.1);
 
    void translate(const double *r);
    void move_to_com();
    void transform_to_principal_axes(int trans_frame=1);
    void transform_to_symmetry_frame();
    void print_xyz(std::ostream& =ExEnv::out0(), const char *title =0) const;
 
    void principal_moments_of_inertia(double *evals, double **evecs=0) const;
 
    int nunique() const { return nuniq_; }
    int unique(int iuniq) const { return equiv_[iuniq][0]; }
    int nequivalent(int iuniq) const { return nequiv_[iuniq]; }
    int equivalent(int iuniq, int j) const { return equiv_[iuniq][j]; }
    int atom_to_unique(int iatom) const { return atom_to_uniq_[iatom]; }
    int atom_to_unique_offset(int iatom) const;
 
    int n_core_electrons();
 
    int max_z();
 
    Ref<AtomInfo> atominfo() const { return atominfo_; }
 
    std::string atom_name(int iatom) const;
 
    std::string atom_symbol(int iatom) const;
 
    void set_include_q(bool iq) { include_q_ = iq; }
    bool include_q() const { return include_q_; }
 
    void set_include_qq(bool iqq) { include_qq_ = iqq; }
    bool include_qq() const { return include_qq_; }
 
    int n_q_atom() const { return q_atoms_.size(); }
    int q_atom(int i) const { return q_atoms_[i]; }
 
    int n_non_q_atom() const { return non_q_atoms_.size(); }
    int non_q_atom(int i) const { return non_q_atoms_[i]; }
 
    bool any_atom_has_charge() const;
    bool any_atom_has_fragment() const;
    bool any_atom_has_label() const;
 
    SCVector3 ref_origin() const { return ref_origin_; }
 
    void save_data_state(StateOut&);
 
#ifdef MPQC_NEW_FEATURES
    virtual boost::property_tree::ptree& write_xml(boost::property_tree::ptree& parent, const XMLWriter& writer);
#endif
 
  private:
    void read_xyz(const char *filename);
#ifdef HAVE_OPENBABEL2
    void read_openbabel2(const char *filename);
#endif // HAVE_OPENBABEL2
};
 
bool operator==(const Molecule& mol1, const Molecule& mol2);
 
// end of addtogroup ChemistryMolecule
 
} // namespace sc
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: