---

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>

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_;

    void clear();
  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::out()) const;
    virtual void print_parsedkeyval(std::ostream& =ExEnv::out(),
                                    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;

    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::out(), 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_; }

    void save_data_state(StateOut&);
};


#endif

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End:

---