Index: ThirdParty/mpqc_open/src/lib/chemistry/molecule/simple.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/molecule/simple.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/molecule/simple.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,253 @@
+
+/* simple.cc -- implementation of the simple internal coordinate classes
+ *
+ *      THIS SOFTWARE FITS THE DESCRIPTION IN THE U.S. COPYRIGHT ACT OF A
+ *      "UNITED STATES GOVERNMENT WORK".  IT WAS WRITTEN AS A PART OF THE
+ *      AUTHOR'S OFFICIAL DUTIES AS A GOVERNMENT EMPLOYEE.  THIS MEANS IT
+ *      CANNOT BE COPYRIGHTED.  THIS SOFTWARE IS FREELY AVAILABLE TO THE
+ *      PUBLIC FOR USE WITHOUT A COPYRIGHT NOTICE, AND THERE ARE NO
+ *      RESTRICTIONS ON ITS USE, NOW OR SUBSEQUENTLY.
+ *
+ *  Author:
+ *      E. T. Seidl
+ *      Bldg. 12A, Rm. 2033
+ *      Computer Systems Laboratory
+ *      Division of Computer Research and Technology
+ *      National Institutes of Health
+ *      Bethesda, Maryland 20892
+ *      Internet: seidl@alw.nih.gov
+ *      February, 1993
+ */
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <string.h>
+#include <math.h>
+
+#if defined(SGI) && !defined(__GNUC__)
+#include <bstring.h>
+#endif
+
+#include <util/class/scexception.h>
+#include <util/misc/formio.h>
+#include <util/state/stateio.h>
+#include <chemistry/molecule/simple.h>
+#include <chemistry/molecule/localdef.h>
+
+#include <util/container/bitarray.h>
+
+using namespace std;
+using namespace sc;
+
+//////////////////////////////////////////////////////////////////////
+
+static ClassDesc SimpleCo_cd(
+  typeid(SimpleCo),"SimpleCo",1,"public IntCoor",
+  0, 0, 0);
+
+SimpleCo::SimpleCo():
+  natoms_(0),
+  atoms(0)
+{
+}
+
+SimpleCo::SimpleCo(int na, const char *re) :
+  IntCoor(re),
+  natoms_(na), atoms(0)
+{
+  atoms=new int[na]; memset(atoms,'\0',sizeof(int)*na);
+}
+
+SimpleCo::SimpleCo(const Ref<KeyVal>&kv,int na) :
+  IntCoor(kv),
+  natoms_(na), atoms(0)
+{
+  atoms=new int[na];
+  memset(atoms,'\0',sizeof(int)*na);
+
+  if (kv->count() == 0) {
+      int i;
+      for (i=0; i<na; i++) {
+          atoms[i] = kv->intvalue("atoms",i);
+          if (kv->error() != KeyVal::OK) break;
+        }
+      if (i == 0) {
+          // couldn't find any atoms so look for a molecule and atom labels
+          Ref<Molecule> mol; mol << kv->describedclassvalue("molecule");
+          if (mol.nonnull()) {
+              for (i=0; i<na; i++) {
+                  char *label = kv->pcharvalue("atom_labels", i);
+                  if (kv->error() != KeyVal::OK) break;
+                  atoms[i] = mol->atom_label_to_index(label) + 1;
+                  delete[] label;
+                  if (atoms[i] == 0) break;
+                }
+            }
+        }
+      if (i != na) {
+          InputError ex("KeyVal CTOR: missing one of the atoms "
+                        "or atom_labels (requires a molecule too) "
+                        "or an atom label was invalid",
+                        __FILE__, __LINE__, 0, 0, class_desc());
+          try {
+              kv->errortrace(ex.elaborate());
+            }
+          catch (...) {}
+          throw ex;
+        }
+    }
+  else {
+      // This is a shorthand form for the input that doesn't allow
+      // the specification of a value.
+      if (label_) delete[] label_;
+      label_=kv->pcharvalue(0);
+      for (int i=0; i<na; i++) {
+          atoms[i]=kv->intvalue(i+1);
+          if (kv->error() != KeyVal::OK) {
+              InputError ex("KeyVal CTOR: missing an atom",
+                            __FILE__, __LINE__, 0, 0, class_desc());
+              try {
+                  kv->errortrace(ex.elaborate());
+                }
+              catch (...) {}
+              throw ex;
+            }
+        }
+    }
+}
+
+SimpleCo::~SimpleCo()
+{
+  if(atoms) delete[] atoms; atoms=0;
+  natoms_=0;
+}
+
+void
+SimpleCo::save_data_state(StateOut& s)
+{
+  IntCoor::save_data_state(s);
+  s.put(natoms_);
+  s.put(atoms,natoms_);
+}
+
+SimpleCo::SimpleCo(StateIn& si):
+  IntCoor(si)
+{
+  si.get(natoms_);
+  si.get(atoms);
+}
+
+int
+SimpleCo::natoms() const
+{
+  return natoms_;
+}
+
+int
+SimpleCo::operator[](int i) const
+{
+  return atoms[i];
+}
+
+int
+SimpleCo::operator!=(SimpleCo&u)
+{
+  return !(*this==u);
+}
+
+int
+SimpleCo::operator==(SimpleCo& sc)
+{
+  if(label_ && !sc.label_ || !label_ && sc.label_) return 0;
+  if(label_ && strcmp(label_,sc.label_)) return 0;
+
+  if(atoms && !sc.atoms || !atoms && sc.atoms) return 0;
+  if(atoms)
+    for(int i=0; i < natoms_; i++) if (atoms[i]!=sc.atoms[i]) return 0;
+
+  return 1;
+}
+
+double
+SimpleCo::force_constant(Ref<Molecule>&mol)
+{
+  return calc_force_con(*mol);
+}
+
+// this updates the values before it computes the bmatrix,
+// which is not quite what I wanted--but close enough
+void
+SimpleCo::bmat(const Ref<Molecule>&mol,RefSCVector&bmat,double coef)
+{
+  int i;
+  int n = bmat.dim().n();
+
+  double* v = new double[n];
+  for (i=0; i<n; i++) v[i] = bmat(i);
+
+  calc_intco(*mol,v,coef);
+
+  for (i=0; i<n; i++) {
+      bmat(i) = v[i];
+    }
+  
+  delete[] v;
+}
+
+void
+SimpleCo::update_value(const Ref<Molecule>&mol)
+{
+  calc_intco(*mol);
+}
+
+void
+SimpleCo::print_details(const Ref<Molecule> &mol, ostream& os) const
+{
+  os << indent
+     << scprintf("%-5s %7s %11.5f", ctype(), (label()?label():""),
+                 preferred_value());
+  
+  int i;
+  for (i=0; i<natoms(); i++)
+      os << scprintf(" %4d", atoms[i]);
+
+  if (mol.nonnull()) {
+      char *separator = " ";
+      os << "  ";
+      for (i=0; i<(4-natoms()); i++) {
+          os << "   ";
+        }
+      for (i=0; i<natoms(); i++) {
+          os << separator << mol->atom_symbol(atoms[i]-1);
+          separator = "-";
+        }
+    }
+
+  os << endl;
+  
+}
+
+// this doesn't catch all cases, it would be best for each subclass
+// to override this
+int
+SimpleCo::equivalent(Ref<IntCoor>&c)
+{
+  if (class_desc() != c->class_desc()) {
+      return 0;
+    }
+  SimpleCo* sc = dynamic_cast<SimpleCo*>(c.pointer());
+  if (natoms_ != sc->natoms_) return 0; // this should never be the case
+  for (int i=0; i<natoms_; i++) {
+      if (atoms[i] != sc->atoms[i]) return 0;
+    }
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
