Index: ThirdParty/mpqc_open/src/lib/chemistry/molecule/molecule.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/molecule/molecule.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/molecule/molecule.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,1433 @@
+//
+// molecule.cc
+//
+// 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.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+
+#include <util/class/scexception.h>
+#include <util/misc/formio.h>
+#include <util/state/stateio.h>
+#include <chemistry/molecule/molecule.h>
+#include <chemistry/molecule/formula.h>
+#include <chemistry/molecule/localdef.h>
+#include <math/scmat/cmatrix.h>
+
+using namespace std;
+using namespace sc;
+
+//////////////////////////////////////////////////////////////////////////
+// Molecule
+
+static ClassDesc Molecule_cd(
+  typeid(Molecule),"Molecule",6,"public SavableState",
+  create<Molecule>, create<Molecule>, create<Molecule>);
+
+Molecule::Molecule():
+  natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
+{
+  pg_ = new PointGroup;
+  atominfo_ = new AtomInfo();
+  geometry_units_ = new Units("bohr");
+  nuniq_ = 0;
+  equiv_ = 0;
+  nequiv_ = 0;
+  atom_to_uniq_ = 0;
+  q_Z_ = atominfo_->string_to_Z("Q");
+  include_q_ = false;
+  include_qq_ = false;
+  init_symmetry_info();
+}
+
+Molecule::Molecule(const Molecule& mol):
+ natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
+{
+  nuniq_ = 0;
+  equiv_ = 0;
+  nequiv_ = 0;
+  atom_to_uniq_ = 0;
+  *this=mol;
+}
+
+Molecule::~Molecule()
+{
+  clear();
+}
+
+void
+Molecule::clear()
+{
+  if (r_) {
+      delete[] r_[0];
+      delete[] r_;
+      r_ = 0;
+    }
+  if (labels_) {
+      for (int i=0; i<natoms_; i++) {
+          delete[] labels_[i];
+        }
+      delete[] labels_;
+      labels_ = 0;
+    }
+  delete[] charges_;
+  charges_ = 0;
+  delete[] mass_;
+  mass_ = 0;
+  delete[] Z_;
+  Z_ = 0;
+
+  clear_symmetry_info();
+}
+
+void
+Molecule::throw_if_atom_duplicated(int begin, double tol)
+{
+  for (int i=begin; i<natoms_; i++) {
+      SCVector3 ri(r_[i]);
+      for (int j=0; j<i; j++) {
+          SCVector3 rj(r_[j]);
+          if (ri.dist(rj) < tol) {
+              throw InputError("duplicated atom coordinate",
+                               __FILE__, __LINE__, 0, 0, class_desc());
+            }
+        }
+    }
+}
+
+Molecule::Molecule(const Ref<KeyVal>&input):
+ natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
+{
+  nuniq_ = 0;
+  equiv_ = 0;
+  nequiv_ = 0;
+  atom_to_uniq_ = 0;
+
+  KeyValValueboolean kvfalse(0);
+  include_q_ = input->booleanvalue("include_q",kvfalse);
+  include_qq_ = input->booleanvalue("include_qq",kvfalse);
+
+  atominfo_ << input->describedclassvalue("atominfo");
+  if (atominfo_.null()) atominfo_ = new AtomInfo;
+  q_Z_ = atominfo_->string_to_Z("Q");
+  if (input->exists("pdb_file")) {
+      geometry_units_ = new Units("angstrom");
+      char* filename = input->pcharvalue("pdb_file");
+      read_pdb(filename);
+      delete[] filename;
+    }
+  else {
+      // check for old style units input first
+      if (input->booleanvalue("angstrom")
+          ||input->booleanvalue("angstroms")
+          ||input->booleanvalue("aangstrom")
+          ||input->booleanvalue("aangstroms")) {
+          geometry_units_ = new Units("angstrom");
+        }
+      // check for new style units input
+      else {
+          geometry_units_ = new Units(input->pcharvalue("unit"),
+                                      Units::Steal);
+        }
+        
+      double conv = geometry_units_->to_atomic_units();
+
+      // get the number of atoms and make sure that the geometry and the
+      // atoms array have the same number of atoms.
+      // right now we read in the unique atoms...then we will symmetrize.
+      // the length of atoms must still equal the length of geometry, but
+      // we'll try to set up atom_labels such that different lengths are
+      // possible
+      int natom = input->count("geometry");
+      if (natom != input->count("atoms")) {
+          throw InputError("size of \"geometry\" != size of \"atoms\"",
+                           __FILE__, __LINE__, 0, 0, class_desc());
+        }
+
+      int i;
+      for (i=0; i<natom; i++) {
+          char *name, *label;
+          int ghost = input->booleanvalue("ghost",i);
+          double charge = input->doublevalue("charge",i);
+          int have_charge = input->error() == KeyVal::OK;
+          if (ghost) {
+              have_charge = 1;
+              charge = 0.0;
+            }
+          add_atom(atominfo_->string_to_Z(name = input->pcharvalue("atoms",i)),
+                   input->doublevalue("geometry",i,0)*conv,
+                   input->doublevalue("geometry",i,1)*conv,
+                   input->doublevalue("geometry",i,2)*conv,
+                   label = input->pcharvalue("atom_labels",i),
+                   input->doublevalue("mass",i),
+                   have_charge, charge
+              );
+          delete[] name;
+          delete[] label;
+        }
+    }
+
+  char *symmetry = input->pcharvalue("symmetry");
+  double symtol = input->doublevalue("symmetry_tolerance",
+                                     KeyValValuedouble(1.0e-4));
+  nuniq_ = 0;
+  equiv_ = 0;
+  nequiv_ = 0;
+  atom_to_uniq_ = 0;
+  if (symmetry && !strcmp(symmetry,"auto")) {
+      set_point_group(highest_point_group(symtol), symtol*10.0);
+    }
+  else {
+      pg_ = new PointGroup(input);
+
+      // translate to the origin of the symmetry frame
+      double r[3];
+      for (int i=0; i<3; i++) {
+          r[i] = -pg_->origin()[i];
+          pg_->origin()[i] = 0;
+        }
+      translate(r);
+
+      if (input->booleanvalue("redundant_atoms")) {
+          init_symmetry_info();
+          cleanup_molecule(symtol);
+        }
+      else {
+          symmetrize();
+          // In case we were given redundant atoms, clean up
+          // the geometry so the symmetry is exact.
+          cleanup_molecule(symtol);
+        }
+    }
+  delete[] symmetry;
+}
+
+Molecule&
+Molecule::operator=(const Molecule& mol)
+{
+  clear();
+
+  pg_ = new PointGroup(*(mol.pg_.pointer()));
+  atominfo_ = mol.atominfo_;
+  geometry_units_ = new Units(mol.geometry_units_->string_rep());
+
+  q_Z_ = mol.q_Z_;
+  include_q_ = mol.include_q_;
+  include_qq_ = mol.include_qq_;
+  q_atoms_ = mol.q_atoms_;
+  non_q_atoms_ = mol.non_q_atoms_;
+
+  natoms_ = mol.natoms_;
+
+  if (natoms_) {
+      if (mol.mass_) {
+          mass_ = new double[natoms_];
+          memcpy(mass_,mol.mass_,natoms_*sizeof(double));
+        }
+      if (mol.charges_) {
+          charges_ = new double[natoms_];
+          memcpy(charges_,mol.charges_,natoms_*sizeof(double));
+        }
+      if (mol.labels_) {
+          labels_ = new char *[natoms_];
+          for (int i=0; i<natoms_; i++) {
+              if (mol.labels_[i]) {
+                  labels_[i] = strcpy(new char[strlen(mol.labels_[i])+1],
+                                      mol.labels_[i]);
+                }
+              else labels_[i] = 0;
+            }
+        }
+      r_ = new double*[natoms_];
+      r_[0] = new double[natoms_*3];
+      for (int i=0; i<natoms_; i++) {
+          r_[i] = &(r_[0][i*3]);
+        }
+      memcpy(r_[0], mol.r_[0], natoms_*3*sizeof(double));
+      Z_ = new int[natoms_];
+      memcpy(Z_, mol.Z_, natoms_*sizeof(int));
+    }
+
+  init_symmetry_info();
+
+  return *this;
+}
+
+void
+Molecule::add_atom(int Z,double x,double y,double z,
+                   const char *label,double mass,
+                   int have_charge, double charge)
+{
+  int i;
+
+  // allocate new arrays
+  int *newZ = new int[natoms_+1];
+  double **newr = new double*[natoms_+1];
+  double *newr0 = new double[(natoms_+1)*3];
+  char **newlabels = 0;
+  if (label || labels_) {
+      newlabels = new char*[natoms_+1];
+    }
+  double *newcharges = 0;
+  if (have_charge || charges_) {
+      newcharges = new double[natoms_+1];
+    }
+  double *newmass = 0;
+  if (mass_ || mass != 0.0) {
+      newmass = new double[natoms_+1];
+    }
+
+  // setup the r_ pointers
+  for (i=0; i<=natoms_; i++) {
+      newr[i] = &(newr0[i*3]);
+    }
+
+  // copy old data to new arrays
+  if (natoms_) {
+      memcpy(newZ,Z_,sizeof(int)*natoms_);
+      memcpy(newr0,r_[0],sizeof(double)*natoms_*3);
+      if (labels_) {
+          memcpy(newlabels,labels_,sizeof(char*)*natoms_);
+        }
+      else if (newlabels) {
+          memset(newlabels,0,sizeof(char*)*natoms_);
+        }
+      if (charges_) {
+          memcpy(newcharges,charges_,sizeof(double)*natoms_);
+        }
+      else if (newcharges) {
+          for (i=0; i<natoms_; i++) newcharges[i] = Z_[i];
+        }
+      if (mass_) {
+          memcpy(newmass,mass_,sizeof(double)*natoms_);
+        }
+      else if (newmass) {
+          memset(newmass,0,sizeof(double)*natoms_);
+        }
+    }
+
+  // delete old data
+  delete[] Z_;
+  if (r_) {
+      delete[] r_[0];
+      delete[] r_;
+    }
+  delete[] labels_;
+  delete[] charges_;
+  delete[] mass_;
+
+  // setup new pointers
+  Z_ = newZ;
+  r_ = newr;
+  labels_ = newlabels;
+  charges_ = newcharges;
+  mass_ = newmass;
+
+  // copy info for this atom into arrays
+  Z_[natoms_] = Z;
+  r_[natoms_][0] = x;
+  r_[natoms_][1] = y;
+  r_[natoms_][2] = z;
+  if (mass_) mass_[natoms_] = mass;
+  if (label) {
+      labels_[natoms_] = strcpy(new char[strlen(label)+1],label);
+    }
+  else if (labels_) {
+      labels_[natoms_] = 0;
+    }
+  if (have_charge) {
+      charges_[natoms_] = charge;
+    }
+  else if (charges_) {
+      charges_[natoms_] = Z;
+    }
+
+  if (Z == q_Z_) {
+      q_atoms_.push_back(natoms_);
+    }
+  else {
+      non_q_atoms_.push_back(natoms_);
+    }
+
+  natoms_++;
+
+  throw_if_atom_duplicated(natoms_-1);
+}
+
+void
+Molecule::print_parsedkeyval(ostream& os,
+                             int print_pg,
+                             int print_unit,
+                             int number_atoms) const
+{
+  int i;
+
+  double conv = geometry_units_->from_atomic_units();
+
+  if (print_pg) pg_->print(os);
+  if (print_unit && geometry_units_->string_rep()) {
+      os << indent
+         << "unit = \"" << geometry_units_->string_rep() << "\""
+         << endl;
+    }
+  os << indent << "{";
+  if (number_atoms) os << scprintf("%3s", "n");
+  os << scprintf(" %5s", "atoms");
+  if (labels_) os << scprintf(" %11s", "atom_labels");
+  int int_charges = 1;
+  if (charges_) {
+      for (i=0;i<natom();i++) if (charges_[i]!=(int)charges_[i]) int_charges=0;
+      if (int_charges) {
+          os << scprintf(" %7s", "charge");
+        }
+      else {
+          os << scprintf(" %17s", "charge");
+        }
+    }
+  os << scprintf("  %16s", "")
+     << scprintf(" %16s", "geometry   ")
+     << scprintf(" %16s ", "");
+  os << "}={" << endl;
+  for (i=0; i<natom(); i++) {
+      os << indent;
+      if (number_atoms) os << scprintf(" %3d", i+1);
+      std::string symbol(atom_symbol(i));
+      os << scprintf(" %5s", symbol.c_str());
+      if (labels_) {
+          const char *lab = labels_[i];
+          if (lab == 0) lab = "";
+          char  *qlab = new char[strlen(lab)+3];
+          strcpy(qlab,"\"");
+          strcat(qlab,lab);
+          strcat(qlab,"\"");
+          os << scprintf(" %11s",qlab);
+          delete[] qlab;
+        }
+      if (charges_) {
+          if (int_charges) os << scprintf(" %7.4f", charges_[i]);
+          else os << scprintf(" %17.15f", charges_[i]);
+        }
+      os << scprintf(" [% 16.10f", conv * r(i,0))
+         << scprintf(" % 16.10f", conv * r(i,1))
+         << scprintf(" % 16.10f]", conv * r(i,2))
+         << endl;
+    }
+  os << indent << "}" << endl;
+}
+
+void
+Molecule::print(ostream& os) const
+{
+  int i;
+
+  MolecularFormula *mf = new MolecularFormula(this);
+  os << indent
+     << "Molecular formula: " << mf->formula() << endl;
+  delete mf;
+
+  os << indent << "molecule<Molecule>: (" << endl;
+  os << incindent;
+  print_parsedkeyval(os);
+  os << decindent;
+  os << indent << ")" << endl;
+
+  os << indent << "Atomic Masses:" << endl;
+  for (i=0; i<natom(); i+=5) {
+      os << indent;
+      for (int j=i; j<i+5 && j<natom(); j++) {
+          os << scprintf(" %10.5f", mass(j));
+        }
+      os << endl;
+    }
+}
+
+int
+Molecule::atom_label_to_index(const char *l) const
+{
+  int i;
+  for (i=0; i<natom(); i++) {
+      if (label(i) && !strcmp(l,label(i))) return i;
+    }
+  return -1;
+}
+
+double*
+Molecule::charges() const
+{
+  double*result = new double[natoms_];
+  if (charges_) {
+      memcpy(result, charges_, sizeof(double)*natom());
+    }
+  else {
+      for (int i=0; i<natom(); i++) result[i] = Z_[i];
+    }
+  return result;
+}
+
+double
+Molecule::charge(int iatom) const
+{
+  if (charges_) return charges_[iatom];
+  return Z_[iatom];
+}
+
+double
+Molecule::nuclear_charge() const
+{
+  double c = 0.0;
+  if (include_q_) {
+    for (int i=0; i<natom(); i++) {
+      c += charge(i);
+    }
+  }
+  else {
+    for (int ii=0; ii<non_q_atoms_.size(); ii++) {
+      c += charge(non_q_atoms_[ii]);
+    }
+  }
+  return c;
+}
+
+void Molecule::save_data_state(StateOut& so)
+{
+  so.put(include_q_);
+  so.put(include_qq_);
+  so.put(natoms_);
+  SavableState::save_state(pg_.pointer(),so);
+  SavableState::save_state(geometry_units_.pointer(),so);
+  SavableState::save_state(atominfo_.pointer(),so);
+  if (natoms_) {
+      so.put(Z_, natoms_);
+      so.put_array_double(r_[0], natoms_*3);
+      so.put(charges_,natoms_);
+    }
+  if (mass_) {
+      so.put(1);
+      so.put_array_double(mass_, natoms_);
+    }
+  else {
+      so.put(0);
+    }
+  if (labels_){
+      so.put(1);
+      for (int i=0; i<natoms_; i++) {
+          so.putstring(labels_[i]);
+        }
+    }
+  else {
+      so.put(0);
+    }
+}
+
+Molecule::Molecule(StateIn& si):
+  SavableState(si),
+  natoms_(0), r_(0), Z_(0), mass_(0), labels_(0)
+{
+  if (si.version(::class_desc<Molecule>()) < 4) {
+      throw FileOperationFailed("cannot restore from old molecules",
+                                __FILE__, __LINE__, 0,
+                                FileOperationFailed::Corrupt,
+                                class_desc());
+    }
+  if (si.version(::class_desc<Molecule>()) < 6) {
+    include_q_ = false;
+    include_qq_ = false;
+    }
+  else {
+    si.get(include_q_);
+    si.get(include_qq_);
+  }
+  si.get(natoms_);
+  pg_ << SavableState::restore_state(si);
+  geometry_units_ << SavableState::restore_state(si);
+  atominfo_ << SavableState::restore_state(si);
+  if (natoms_) {
+      si.get(Z_);
+      r_ = new double*[natoms_];
+      r_[0] = new double[natoms_*3];
+      si.get_array_double(r_[0],natoms_*3);
+      for (int i=1; i<natoms_; i++) {
+          r_[i] = &(r_[0][i*3]);
+        }
+      if (si.version(::class_desc<Molecule>()) > 4) {
+          si.get(charges_);
+        }
+      else {
+          charges_ = 0;
+        }
+    }
+  int test;
+  si.get(test);
+  if (test) {
+      mass_ = new double[natoms_];
+      si.get_array_double(mass_, natoms_);
+    }
+  si.get(test);
+  if (test){
+      labels_ = new char*[natoms_];
+      for (int i=0; i<natoms_; i++) {
+          si.getstring(labels_[i]);
+        }
+    }
+
+  for (int i=0; i<natoms_; i++) {
+      if (Z_[i] == q_Z_) {
+          q_atoms_.push_back(i);
+        }
+      else {
+          non_q_atoms_.push_back(i);
+        }
+    }
+
+  nuniq_ = 0;
+  equiv_ = 0;
+  nequiv_ = 0;
+  atom_to_uniq_ = 0;
+  init_symmetry_info();
+}
+
+int
+Molecule::atom_to_unique_offset(int iatom) const
+{
+  int iuniq = atom_to_uniq_[iatom];
+  int nequiv = nequiv_[iuniq];
+  for (int i=0; i<nequiv; i++) {
+      if (equiv_[iuniq][i] == iatom) return i;
+    }
+  ExEnv::errn() << "Molecule::atom_to_unique_offset: internal error"
+               << endl;
+  return -1;
+}
+
+void
+Molecule::set_point_group(const Ref<PointGroup>&ppg, double tol)
+{
+  ExEnv::out0() << indent
+       << "Molecule: setting point group to " << ppg->symbol()
+       << endl;
+  pg_ = new PointGroup(*ppg.pointer());
+
+  double r[3];
+  for (int i=0; i<3; i++) {
+      r[i] = -pg_->origin()[i];
+      pg_->origin()[i] = 0;
+    }
+  translate(r);
+
+  clear_symmetry_info();
+  init_symmetry_info();
+
+  cleanup_molecule(tol);
+}
+
+Ref<PointGroup>
+Molecule::point_group() const
+{
+  return pg_;
+}
+
+SCVector3
+Molecule::center_of_mass() const
+{
+  SCVector3 ret;
+  double M;
+
+  ret = 0.0;
+  M = 0.0;
+
+  for (int i=0; i < natom(); i++) {
+    double m = mass(i);
+    ret += m * SCVector3(r(i));
+    M += m;
+  }
+
+  ret *= 1.0/M;
+
+  return ret;
+}
+
+double
+Molecule::nuclear_repulsion_energy()
+{
+  double e=0.0;
+
+  // non_q non_q terms
+  for (int ii=1; ii < non_q_atoms_.size(); ii++) {
+      int i = non_q_atoms_[ii];
+      SCVector3 ai(r(i));
+      double Zi = charge(i);
+    
+      for (int jj=0; jj < ii; jj++) {
+          int j = non_q_atoms_[jj];
+          SCVector3 aj(r(j));
+          e += Zi * charge(j) / ai.dist(aj);
+        }
+    }
+
+  // non_q q terms
+  for (int ii=0; ii < q_atoms_.size(); ii++) {
+      int i = q_atoms_[ii];
+      SCVector3 ai(r(i));
+      double Zi = charge(i);
+    
+      for (int jj=0; jj < non_q_atoms_.size(); jj++) {
+          int j = non_q_atoms_[jj];
+          SCVector3 aj(r(j));
+          e += Zi * charge(j) / ai.dist(aj);
+        }
+    }
+
+  if (include_qq_) {
+      // q q terms
+      for (int ii=1; ii < q_atoms_.size(); ii++) {
+          int i = q_atoms_[ii];
+          SCVector3 ai(r(i));
+          double Zi = charge(i);
+    
+          for (int jj=0; jj < ii; jj++) {
+              int j = q_atoms_[jj];
+              SCVector3 aj(r(j));
+              e += Zi * charge(j) / ai.dist(aj);
+            }
+        }
+    }
+
+  return e;
+}
+
+void
+Molecule::nuclear_repulsion_1der(int center, double xyz[3])
+{
+  int i,j,k;
+  double rd[3],r2;
+  double factor;
+
+  xyz[0] = 0.0;
+  xyz[1] = 0.0;
+  xyz[2] = 0.0;
+
+  SCVector3 r_center(r(center));
+  double Z_center = charge(center);
+  bool center_is_Q = (atom_symbol(center) == "Q");
+
+  // this handles center = Q or non_Q and atom = non_Q
+  for (int ii=0; ii < non_q_atoms_.size(); ii++) {
+    int i = non_q_atoms_[ii];
+    if (i == center) continue;
+    SCVector3 r_i(r(i));
+
+    r2 = 0.0;
+    for (k=0; k < 3; k++) {
+      rd[k] = r_center[k] - r_i[k];
+      r2 += rd[k]*rd[k];
+    }
+    factor = - Z_center * charge(i) * pow(r2,-1.5);
+    for (k=0; k<3; k++) {
+      xyz[k] += factor * rd[k];
+    }
+  }
+
+  // this handles center = Q or non_Q and atom = Q
+  for (int ii=0; ii < q_atoms_.size(); ii++) {
+    int i = q_atoms_[ii];
+    if (i == center || (!include_qq_ && center_is_Q)) continue;
+    SCVector3 r_i(r(i));
+
+    r2 = 0.0;
+    for (k=0; k < 3; k++) {
+      rd[k] = r_center[k] - r_i[k];
+      r2 += rd[k]*rd[k];
+    }
+    factor = - Z_center * charge(i) * pow(r2,-1.5);
+    for (k=0; k<3; k++) {
+      xyz[k] += factor * rd[k];
+    }
+  }
+}
+
+void
+Molecule::nuclear_charge_efield(const double *charges,
+                                const double *position, double *efield)
+{
+  double tmp;
+  double rd[3];
+
+  for (int i=0; i<3; i++) efield[i] = 0.0;
+
+  if (include_q_) {
+    for (int i=0; i<natoms_; i++) {
+      SCVector3 a(r(i));
+      tmp = 0.0;
+      for (int j=0; j<3; j++) {
+        rd[j] = position[j] - a[j];
+        tmp += rd[j]*rd[j];
+      }
+      tmp = charges[i]/(tmp*sqrt(tmp));
+      for (int j=0; j<3; j++) {
+        efield[j] +=  rd[j] * tmp;
+      }
+    }
+  }
+  else {
+    for (int ii=0; ii<non_q_atoms_.size(); ii++) {
+      int i = non_q_atoms_[ii];
+      SCVector3 a(r(i));
+      tmp = 0.0;
+      for (int j=0; j<3; j++) {
+        rd[j] = position[j] - a[j];
+        tmp += rd[j]*rd[j];
+      }
+      tmp = charges[i]/(tmp*sqrt(tmp));
+      for (int j=0; j<3; j++) {
+        efield[j] +=  rd[j] * tmp;
+      }
+    }
+  }
+}
+
+void
+Molecule::nuclear_efield(const double *position, double *efield)
+{
+  double tmp;
+  double rd[3];
+
+  for (int i=0; i<3; i++) efield[i] = 0.0;
+
+  if (include_q_) {
+    for (int i=0; i<natoms_; i++) {
+      SCVector3 a(r(i));
+      tmp = 0.0;
+      for (int j=0; j<3; j++) {
+        rd[j] = position[j] - a[j];
+        tmp += rd[j]*rd[j];
+      }
+      tmp = charge(i)/(tmp*sqrt(tmp));
+      for (int j=0; j<3; j++) {
+        efield[j] +=  rd[j] * tmp;
+      }
+    }
+  }
+  else {
+    for (int ii=0; ii<non_q_atoms_.size(); ii++) {
+      int i = non_q_atoms_[ii];
+      SCVector3 a(r(i));
+      tmp = 0.0;
+      for (int j=0; j<3; j++) {
+        rd[j] = position[j] - a[j];
+        tmp += rd[j]*rd[j];
+      }
+      tmp = charge(i)/(tmp*sqrt(tmp));
+      for (int j=0; j<3; j++) {
+        efield[j] +=  rd[j] * tmp;
+      }
+    }
+  }
+}
+
+int
+Molecule::atom_at_position(double *v, double tol) const
+{
+  SCVector3 p(v);
+  for (int i=0; i < natom(); i++) {
+      SCVector3 ai(r(i));
+      if (p.dist(ai) < tol) return i;
+    }
+  return -1;
+}
+
+void
+Molecule::symmetrize(const Ref<PointGroup> &pg, double tol)
+{
+  pg_ = new PointGroup(pg);
+
+  // translate to the origin of the symmetry frame
+  double r[3];
+  for (int i=0; i<3; i++) {
+      r[i] = -pg_->origin()[i];
+      pg_->origin()[i] = 0;
+    }
+  translate(r);
+
+  symmetrize(tol);
+}
+
+// We are given a molecule which may or may not have just the symmetry
+// distinct atoms in it.  We have to go through the existing set of atoms,
+// perform each symmetry operation in the point group on each of them, and
+// then add the new atom if it isn't in the list already
+
+void
+Molecule::symmetrize(double tol)
+{
+  // if molecule is c1, don't do anything
+  if (!strcmp(this->point_group()->symbol(),"c1")) {
+    init_symmetry_info();
+    return;
+    }
+
+  clear_symmetry_info();
+
+  Molecule *newmol = new Molecule(*this);
+  
+  CharacterTable ct = this->point_group()->char_table();
+
+  SCVector3 np;
+  SymmetryOperation so;
+  
+  for (int i=0; i < natom(); i++) {
+    SCVector3 ac(r(i));
+
+    for (int g=0; g < ct.order(); g++) {
+      so = ct.symm_operation(g);
+      for (int ii=0; ii < 3; ii++) {
+        np[ii]=0;
+        for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * ac[jj];
+      }
+
+      int atom = newmol->atom_at_position(np.data(), tol);
+      if (atom < 0) {
+        newmol->add_atom(Z_[i],np[0],np[1],np[2],label(i));
+      }
+      else {
+        if (Z(i) != newmol->Z(atom)
+            || fabs(mass(i)-newmol->mass(atom))>1.0e-10) {
+            throw ToleranceExceeded("symmetrize: atom mismatch",
+                                    __FILE__, __LINE__,
+                                    1.0e-10, fabs(mass(i)-newmol->mass(atom)),
+                                    class_desc());
+        }
+      }
+    }
+  }
+  
+  Ref<Units> saved_units = geometry_units_;
+  *this = *newmol;
+  geometry_units_ = saved_units;
+  delete newmol;
+
+  init_symmetry_info();
+}
+
+void
+Molecule::translate(const double *r)
+{
+  for (int i=0; i < natom(); i++) {
+    r_[i][0] += r[0];
+    r_[i][1] += r[1];
+    r_[i][2] += r[2];
+  }
+}
+
+// move the molecule to the center of mass
+void
+Molecule::move_to_com()
+{
+  SCVector3 com = -center_of_mass();
+  translate(com.data());
+}
+
+// find the 3 principal coordinate axes, and rotate the molecule to be 
+// aligned along them.  also rotate the symmetry frame contained in point_group
+void
+Molecule::transform_to_principal_axes(int trans_frame)
+{
+  // mol_move_to_com(mol);
+
+  double *inert[3], inert_dat[9], *evecs[3], evecs_dat[9];
+  double evals[3];
+
+  int i,j,k;
+  for (i=0; i < 3; i++) {
+    inert[i] = &inert_dat[i*3];
+    evecs[i] = &evecs_dat[i*3];
+  }
+  memset(inert_dat,'\0',sizeof(double)*9);
+  memset(evecs_dat,'\0',sizeof(double)*9);
+
+  for (i=0; i < natom(); i++) {
+    SCVector3 ac(r(i));
+    double m=mass(i);
+    inert[0][0] += m * (ac[1]*ac[1] + ac[2]*ac[2]);
+    inert[1][0] -= m * ac[0]*ac[1];
+    inert[1][1] += m * (ac[0]*ac[0] + ac[2]*ac[2]);
+    inert[2][0] -= m * ac[0]*ac[2];
+    inert[2][1] -= m * ac[1]*ac[2];
+    inert[2][2] += m * (ac[0]*ac[0] + ac[1]*ac[1]);
+  }
+
+  inert[0][1] = inert[1][0];
+  inert[0][2] = inert[2][0];
+  inert[1][2] = inert[2][1];
+
+ // cleanup inert
+  for (i=0; i < 3; i++) {
+    for (int j=0; j <= i; j++) {
+      if (fabs(inert[i][j]) < 1.0e-5) {
+        inert[i][j]=inert[j][i]=0.0;
+      }
+    }
+  }
+
+  cmat_diag(inert, evals, evecs, 3, 1, 1e-14);
+
+ // cleanup evecs
+  for (i=0; i < 3; i++) {
+    for (int j=0; j < 3; j++) {
+      if (fabs(evecs[i][j]) < 1.0e-5) {
+        evecs[i][j]=0.0;
+      }
+    }
+  }
+
+  for (i=0; i<natom(); i++) {
+      double a[3];
+      a[0] = r(i,0); a[1] = r(i,1); a[2] = r(i,2);
+      for (j=0; j<3; j++) {
+          double e = 0.0;
+          for (k=0; k<3; k++) {
+              e += a[k] * evecs[k][j];
+            }
+          r_[i][j] = e;
+        }
+    }
+
+  if (!trans_frame) return;
+  
+  SymmetryOperation tso=point_group()->symm_frame();
+
+  for (i=0; i < 3; i++) {
+    for (int j=0; j < 3; j++) {
+      double t=0;
+      for (int k=0; k < 3; k++) t += tso[k][j]*evecs[k][i];
+      pg_->symm_frame()[i][j] = t;
+    }
+  }
+}
+
+void
+Molecule::transform_to_symmetry_frame()
+{
+  int i,j,k;
+  double t[3][3];
+
+  SymmetryOperation tso=point_group()->symm_frame();
+
+  for (i=0; i<3; i++) {
+      for (j=0; j<3; j++) {
+          t[i][j] = tso[i][j];
+        }
+    }
+
+  for (i=0; i<natom(); i++) {
+      double a[3];
+      a[0] = r(i,0); a[1] = r(i,1); a[2] = r(i,2);
+      for (j=0; j<3; j++) {
+          double e = 0.0;
+          for (k=0; k<3; k++) {
+              e += a[k] * t[k][j];
+            }
+          r_[i][j] = e;
+        }
+    }
+
+  for (i=0; i<3; i++) {
+      for (j=0; j<3; j++) {
+          double e=0;
+          for (k=0; k<3; k++) e += tso[k][j]*t[k][i];
+          pg_->symm_frame()[i][j] = e;
+        }
+    }
+}
+
+// given a molecule, make sure that equivalent centers have coordinates
+// that really map into each other
+
+void
+Molecule::cleanup_molecule(double tol)
+{
+  // if symmetry is c1, do nothing else
+  if (!strcmp(point_group()->symbol(),"c1")) return;
+
+  int i;
+  SCVector3 up,np,ap;
+  SymmetryOperation so;
+  CharacterTable ct = point_group()->char_table();
+
+  // first clean up the unique atoms by replacing each coordinate with the
+  // average of coordinates obtained by applying all symmetry operations to
+  // the original atom, iff the new atom ends up near the original atom
+  for (i=0; i < nunique(); i++) {
+      // up will store the original coordinates of unique atom i
+      up = r(unique(i));
+      // ap will hold the average coordinate (times the number of coordinates)
+      // initialize it to the E result
+      ap = up;
+      int ncoor = 1;
+      // loop through all sym ops except E
+      for (int g=1; g < ct.order(); g++) {
+          so = ct.symm_operation(g);
+          for (int ii=0; ii < 3; ii++) {
+              np[ii]=0;
+              for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * up[jj];
+            }
+          if (np.dist(up) < 0.1) {
+              for (int jj=0; jj < 3; jj++) ap[jj] += np[jj];
+              ncoor++;
+            }
+        }
+      // replace the unique coordinate with the average coordinate
+      r_[unique(i)][0] = ap[0] / ncoor;
+      r_[unique(i)][1] = ap[1] / ncoor;
+      r_[unique(i)][2] = ap[2] / ncoor;
+    }
+
+  // find the atoms equivalent to each unique atom and eliminate
+  // numerical errors that may be in the equivalent atom's coordinates
+
+  // loop through unique atoms
+  for (i=0; i < nunique(); i++) {
+      // up will store the coordinates of unique atom i
+      up = r(unique(i));
+
+      // loop through all sym ops except E
+      for (int g=1; g < ct.order(); g++) {
+          so = ct.symm_operation(g);
+          for (int ii=0; ii < 3; ii++) {
+              np[ii]=0;
+              for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * up[jj];
+            }
+
+          // loop through equivalent atoms
+          int found = 0;
+          for (int j=0; j < natom(); j++) {
+              // see if j is generated from i
+              if (np.dist(SCVector3(r(j))) < tol) {
+                  r_[j][0] = np[0];
+                  r_[j][1] = np[1];
+                  r_[j][2] = np[2];
+                  found = 1;
+                }
+            }
+          if (!found) {
+              SCException ex("cleanup: couldn't find atom",
+                             __FILE__, __LINE__, class_desc());
+              try {
+                  ex.elaborate()
+                      << "couldn't find atom at " << np << endl
+                      << "transforming uniq atom " << i << " at " << up << endl
+                      << "with symmetry op " << g << ":" << endl;
+                  so.print(ex.elaborate());
+                }
+              catch (...) {}
+              throw ex;
+            }
+        }
+    }
+
+}
+
+///////////////////////////////////////////////////////////////////
+// Compute the principal axes and the principal moments of inertia
+///////////////////////////////////////////////////////////////////
+
+void
+Molecule::principal_moments_of_inertia(double *evals, double **evecs) const
+{
+
+  // The principal moments of inertia are computed in amu*angstrom^2
+  // evals: principal moments of inertia
+  // evecs: principal axes (optional argument)
+
+  Ref<Units> units = new Units("angstroms * angstroms");
+  double au_to_angs = units->from_atomic_units();
+
+  double *inert[3];  // inertia tensor
+
+  int i, j;
+  int delete_evecs = 0;
+
+  // (allocate and) initialize evecs, evals, and inert
+  if (!evecs) {
+    evecs = new double*[3];
+    for (i=0; i<3; i++) evecs[i] = new double[3];
+    delete_evecs = 1;
+    }
+  for (i=0; i<3; i++) {
+    inert[i] = new double[3];
+    memset(inert[i],'\0',sizeof(double)*3);
+    memset(evecs[i],'\0',sizeof(double)*3);
+    }
+  memset(evals,'\0',sizeof(double)*3);
+
+  SCVector3 com = center_of_mass();
+
+  // compute inertia tensor
+  SCVector3 ac;
+  for (i=0; i<natom(); i++) {
+    ac = r(i);
+    // compute moments of inertia wrt center of mass
+    for (j=0; j<3; j++) ac(j) -= com(j);
+    double m=au_to_angs*mass(i);
+    inert[0][0] += m * (ac[1]*ac[1] + ac[2]*ac[2]);
+    inert[1][0] -= m * ac[0]*ac[1];
+    inert[1][1] += m * (ac[0]*ac[0] + ac[2]*ac[2]);
+    inert[2][0] -= m * ac[0]*ac[2];
+    inert[2][1] -= m * ac[1]*ac[2];
+    inert[2][2] += m * (ac[0]*ac[0] + ac[1]*ac[1]);
+    }
+  inert[0][1] = inert[1][0];
+  inert[0][2] = inert[2][0];
+  inert[1][2] = inert[2][1];
+
+  cmat_diag(inert, evals, evecs, 3, 1, 1e-14);
+
+  if (delete_evecs) {
+    for (i=0; i<3; i++) delete[] evecs[i];
+    delete[] evecs;
+    }
+  for (i=0; i<3; i++) {
+    delete[] inert[i];
+    }
+}
+
+int
+Molecule::n_core_electrons()
+{
+  int i,n=0;
+  for (i=0; i<natom(); i++) {
+      if (charge(i) == 0.0) continue;
+      int z = Z_[i];
+      if (z > 2) n += 2;
+      if (z > 10) n += 8;
+      if (z > 18) n += 8;
+      if (z > 30) n += 10;
+      if (z > 36) n += 8;
+      if (z > 48) n += 10;
+      if (z > 54) n += 8;
+      if (z > 72) {
+          throw LimitExceeded<int>("n_core_electrons: atomic number too large",
+                                   __FILE__, __LINE__, 72, z, class_desc());
+        }
+    }
+  return n;
+}
+
+int
+Molecule::max_z()
+{
+  int i, maxz=0;
+  for (i=0; i<natom(); i++) {
+      int z = Z_[i];
+      if (z>maxz) maxz = z;
+    }
+  return maxz;
+}
+
+void
+Molecule::read_pdb(const char *filename)
+{
+  clear();
+  ifstream in(filename);
+  Ref<Units> units = new Units("angstrom");
+  while (in.good()) {
+      const int max_line = 80;
+      char line[max_line];
+      in.getline(line,max_line);
+      char *endofline = (char*) memchr(line, 0, max_line);
+      if (endofline) memset(endofline, ' ', &line[max_line-1] - endofline);
+      if (!in.good()) break;
+      if (strncmp(line,"ATOM  ",6) == 0
+          ||strncmp(line,"HETATM",6) == 0) {
+
+          char element[3];
+          strncpy(element,&line[76],2); element[2] = '\0';
+          char name[5];
+          strncpy(name,&line[12],4); name[4] = '\0';
+          if (element[0]==' '&&element[1]==' ') {
+              // no element was given so get the element from the atom name
+              if (name[0]!=' '&&name[3]!=' ') {
+                  // some of the atom label may have been
+                  // pushed over into the element fields
+                  // so check the residue
+                  char resName[4];
+                  strncpy(resName,&line[17],3); resName[3] = '\0';
+                  if (strncmp(line,"ATOM  ",6)==0&&(0
+                      ||strcmp(resName,"ALA")==0||strcmp(resName,"A  ")==0
+                      ||strcmp(resName,"ARG")==0||strcmp(resName,"R  ")==0
+                      ||strcmp(resName,"ASN")==0||strcmp(resName,"N  ")==0
+                      ||strcmp(resName,"ASP")==0||strcmp(resName,"D  ")==0
+                      ||strcmp(resName,"ASX")==0||strcmp(resName,"B  ")==0
+                      ||strcmp(resName,"CYS")==0||strcmp(resName,"C  ")==0
+                      ||strcmp(resName,"GLN")==0||strcmp(resName,"Q  ")==0
+                      ||strcmp(resName,"GLU")==0||strcmp(resName,"E  ")==0
+                      ||strcmp(resName,"GLX")==0||strcmp(resName,"Z  ")==0
+                      ||strcmp(resName,"GLY")==0||strcmp(resName,"G  ")==0
+                      ||strcmp(resName,"HIS")==0||strcmp(resName,"H  ")==0
+                      ||strcmp(resName,"ILE")==0||strcmp(resName,"I  ")==0
+                      ||strcmp(resName,"LEU")==0||strcmp(resName,"L  ")==0
+                      ||strcmp(resName,"LYS")==0||strcmp(resName,"K  ")==0
+                      ||strcmp(resName,"MET")==0||strcmp(resName,"M  ")==0
+                      ||strcmp(resName,"PHE")==0||strcmp(resName,"F  ")==0
+                      ||strcmp(resName,"PRO")==0||strcmp(resName,"P  ")==0
+                      ||strcmp(resName,"SER")==0||strcmp(resName,"S  ")==0
+                      ||strcmp(resName,"THR")==0||strcmp(resName,"T  ")==0
+                      ||strcmp(resName,"TRP")==0||strcmp(resName,"W  ")==0
+                      ||strcmp(resName,"TYR")==0||strcmp(resName,"Y  ")==0
+                      ||strcmp(resName,"VAL")==0||strcmp(resName,"V  ")==0
+                      ||strcmp(resName,"A  ")==0
+                      ||strcmp(resName,"+A ")==0
+                      ||strcmp(resName,"C  ")==0
+                      ||strcmp(resName,"+C ")==0
+                      ||strcmp(resName,"G  ")==0
+                      ||strcmp(resName,"+G ")==0
+                      ||strcmp(resName,"I  ")==0
+                      ||strcmp(resName,"+I ")==0
+                      ||strcmp(resName,"T  ")==0
+                      ||strcmp(resName,"+T ")==0
+                      ||strcmp(resName,"U  ")==0
+                      ||strcmp(resName,"+U ")==0
+                      )) {
+                      // there no two letter elements for these cases
+                      element[0] = name[0];
+                      element[1] = '\0';
+                    }
+                }
+              else {
+                  strncpy(element,name,2); element[2] = '\0';
+                }
+            }
+          if (element[0] == ' ') {
+              element[0] = element[1];
+              element[1] = '\0';
+            }
+
+          int Z = atominfo_->string_to_Z(element);
+
+          char field[9];
+          strncpy(field,&line[30],8); field[8] = '\0';
+          double x = atof(field);
+          strncpy(field,&line[38],8); field[8] = '\0';
+          double y = atof(field);
+          strncpy(field,&line[46],8); field[8] = '\0';
+          double z = atof(field);
+          add_atom(Z,
+                   x*units->to_atomic_units(),
+                   y*units->to_atomic_units(),
+                   z*units->to_atomic_units());
+        }
+      else {
+        // skip to next record
+        }
+    }
+}
+
+void
+Molecule::print_pdb(ostream& os, char *title) const
+{
+  Ref<Units> u = new Units("angstrom");
+  double bohr = u->from_atomic_units();
+
+  if (title)
+      os << scprintf("%-10s%-60s\n","COMPND",title);
+  else
+      os << scprintf("%-10s%-60s\n","COMPND","Title");
+
+  if (title)
+      os << scprintf("REMARK   %s\n", title);
+
+  int i;
+  for (i=0; i < natom(); i++) {
+    char symb[4];
+    std::string symbol(atom_symbol(i));
+    sprintf(symb,"%s1",symbol.c_str());
+
+    os << scprintf(
+        "HETATM%5d  %-3s UNK %5d    %8.3f%8.3f%8.3f  0.00  0.00   0\n",
+        i+1, symb, 0, r(i,0)*bohr, r(i,1)*bohr, r(i,2)*bohr);
+  }
+
+  for (i=0; i < natom(); i++) {
+    double at_rad_i = atominfo_->atomic_radius(Z_[i]);
+    SCVector3 ai(r(i));
+
+    os << scprintf("CONECT%5d",i+1);
+
+    for (int j=0; j < natom(); j++) {
+
+      if (j==i) continue;
+
+      double at_rad_j = atominfo_->atomic_radius(Z_[j]);
+      SCVector3 aj(r(j));
+
+      if (ai.dist(aj) < 1.1*(at_rad_i+at_rad_j))
+          os << scprintf("%5d",j+1);
+    }
+
+    os << endl;
+  }
+
+  os << "END" << endl;
+  os.flush();
+}
+
+double
+Molecule::mass(int atom) const
+{
+  if (!mass_ || mass_[atom] == 0) {
+      return atominfo_->mass(Z_[atom]);
+    }
+  return mass_[atom];
+}
+
+const char *
+Molecule::label(int atom) const
+{
+  if (!labels_) return 0;
+  return labels_[atom];
+}
+
+std::string
+Molecule::atom_name(int iatom) const
+{
+  return atominfo_->name(Z_[iatom]);
+}
+
+std::string
+Molecule::atom_symbol(int iatom) const
+{
+  return atominfo_->symbol(Z_[iatom]);
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
