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
#ifdef __GNUC__
#pragma interface
#endif
#include <stdio.h>
#include <iostream>
#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>
namespace sc {
class Molecule: public SavableState
{
protected:
int natoms_;
Ref<AtomInfo> atominfo_;
Ref<PointGroup> pg_;
Ref<Units> geometry_units_;
double **r_;
int *Z_;
double *charges_;
// 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();
// these are optional
double *mass_;
char **labels_;
// 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 char * = 0, double mass = 0.0,
int have_charge = 0, double charge = 0.0);
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;
int natom() const { return natoms_; }
int Z(int atom) const { return Z_[atom]; }
double &r(int atom, int xyz) { return r_[atom][xyz]; }
const double &r(int atom, int xyz) const { return r_[atom][xyz]; }
double *r(int atom) { return r_[atom]; }
const double *r(int atom) const { return r_[atom]; }
double mass(int atom) const;
const char *label(int atom) const;
int atom_at_position(double *, double tol = 0.05) const;
int atom_label_to_index(const char *label) const;
double *charges() const;
double charge(int iatom) const;
double nuclear_charge() const;
void set_point_group(const Ref<PointGroup>&, double tol=1.0e-7);
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;
SCVector3 center_of_charge() 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_pdb(std::ostream& =ExEnv::out0(), char *title =0) const;
void read_pdb(const char *filename);
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]; }
void save_data_state(StateOut&);
};
}
#endif
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: