Index: ThirdParty/mpqc_open/src/lib/chemistry/molecule/atominfo.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/molecule/atominfo.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/molecule/atominfo.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,611 @@
+//
+// molinfo.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.
+//
+
+#include <stdlib.h>
+#include <util/misc/string.h>
+#include <ctype.h>
+#include <sys/stat.h>
+
+#include <sstream>
+#include <stdexcept>
+
+#include <util/misc/units.h>
+#include <util/misc/autovec.h>
+#include <util/misc/formio.h>
+#include <util/state/stateio.h>
+#include <util/group/message.h>
+#include <chemistry/molecule/atominfo.h>
+
+using namespace std;
+using namespace sc;
+
+////////////////////////////////////////////////////////////////////////
+// AtomInfo
+
+struct AtomInfo::atom
+AtomInfo::elements_[Nelement] = 
+  {{1, "hydrogen",   "H"},
+   {2, "helium",     "He"},
+   {3, "lithium",    "Li"},
+   {4, "beryllium",  "Be"},
+   {5, "boron",      "B"},
+   {6, "carbon",     "C"},
+   {7, "nitrogen",   "N"},
+   {8, "oxygen",     "O"},
+   {9, "fluorine",   "F"},
+   {10, "neon",       "Ne"},
+   {11, "sodium",     "Na"},
+   {12, "magnesium",  "Mg"},
+   {13, "aluminum",   "Al"},
+   {14, "silicon",    "Si"},
+   {15, "phosphorus" ,"P"},
+   {16, "sulfur",     "S"},
+   {17, "chlorine",   "Cl"},
+   {18, "argon",      "Ar"},
+   {19, "potassium",  "K"},
+   {20, "calcium",    "Ca"},
+   {21, "scandium",   "Sc"},
+   {22, "titanium",   "Ti"},
+   {23, "vanadium",   "V"},
+   {24, "chromium",   "Cr"},
+   {25, "manganese",  "Mn"},
+   {26, "iron",       "Fe"},
+   {27, "cobalt",     "Co"},
+   {28, "nickel",     "Ni"},
+   {29, "copper",     "Cu"},
+   {30, "zinc",       "Zn"},
+   {31, "gallium",    "Ga"},
+   {32, "germanium",  "Ge"},
+   {33, "arsenic",    "As"},
+   {34, "selenium",   "Se"},
+   {35, "bromine",    "Br"},
+   {36, "krypton",    "Kr"},
+   {37, "rubidium",     "Rb"},
+   {38, "strontium",    "Sr"},
+   {39, "yttrium",      "Y"},
+   {40, "zirconium",    "Zr"},
+   {41, "niobium",      "Nb"},
+   {42, "molybdenum",   "Mo"},
+   {43, "technetium",   "Tc"},
+   {44, "ruthenium",    "Ru"},
+   {45, "rhodium",      "Rh"},
+   {46, "palladium",    "Pd"},
+   {47, "silver",       "Ag"},
+   {48, "cadminium",    "Cd"},
+   {49, "indium",       "In"},
+   {50, "tin",          "Sn"},
+   {51, "antimony",     "Sb"},
+   {52, "tellurium",    "Te"},
+   {53, "iodine",       "I"},
+   {54, "xenon",        "Xe"},
+   {55, "cesium",       "Cs"},
+   {56, "barium",       "Ba"},
+   {57, "lanthanium",   "La"},
+   {58, "cerium",       "Ce"},
+   {59, "praseodymium", "Pr"},
+   {60, "neodymium",    "Nd"},
+   {61, "promethium",   "Pm"},
+   {62, "samarium",     "Sm"},
+   {63, "europium",     "Eu"},
+   {64, "gadolinium",   "Gd"},
+   {65, "terbium",      "Tb"},
+   {66, "dysprosium",   "Dy"},
+   {67, "holmium",      "Ho"},
+   {68, "erbium",       "Er"},
+   {69, "thulium",      "Tm"},
+   {70, "ytterbium",    "Yb"},
+   {71, "lutetium",     "Lu"},
+   {72, "hafnium",      "Hf"},
+   {73, "tantalum",     "Ta"},
+   {74, "tungsten",     "W"},
+   {75, "rhenium",      "Re"},
+   {76, "osmium",       "Os"},
+   {77, "iridium",      "Ir"},
+   {78, "platinum",     "Pt"},
+   {79, "gold",         "Au"},
+   {80, "mercury",      "Hg"},
+   {81, "thallium",     "Tl"},
+   {82, "lead",         "Pb"},
+   {83, "bismuth",      "Bi"},
+   {84, "polonium",     "Po"},
+   {85, "astatine",     "At"},
+   {86, "radon",        "Rn"},
+   {87, "francium",     "Fr"},
+   {88, "radium",       "Ra"},
+   {89, "actinium",     "Ac"},
+   {90, "thorium",      "Th"},
+   {91, "protactinium", "Pa"},
+   {92, "uranium",      "U"},
+   {93, "neptunium",    "Np"},
+   {94, "plutonium",    "Pu"},
+   {95, "americium",    "Am"},
+   {96, "curium",       "Cm"},
+   {97, "berkelium",    "Bk"},
+   {98, "californium",  "Cf"},
+   {99, "einsteinum",   "Es"},
+   {100, "fermium",      "Fm"},
+   {101, "mendelevium",  "Md"},
+   {102, "nobelium",     "No"},
+   {103, "lawrencium",   "Lr"},
+   {104, "rutherfordium","Rf"},
+   {105, "hahnium",      "Ha"},
+   {106, "seaborgium",   "Sg"},
+   {107, "bohrium",      "Bh"},
+   {108, "hassium",      "Hs"},
+   {109, "meitnerium",   "Mt"},
+   {110, "darmstadtium", "Ds"},
+   {111, "roentgenium",  "Rg"},
+   {112, "ununbium",     "Uub"},
+   {113, "ununtrium",    "Uut"},
+   {114, "ununquadium",  "Uuq"},
+   {115, "ununpentium",  "Uup"},
+   {116, "ununhexium",   "Uuh"},
+   {117, "ununseptium",  "Uus"},
+   {118, "ununoctium",   "Uuo"}
+  };
+
+static ClassDesc AtomInfo_cd(
+  typeid(AtomInfo),"AtomInfo",3,"public SavableState",
+  0, create<AtomInfo>, create<AtomInfo>);
+
+AtomInfo::AtomInfo()
+{
+  initialize_names();
+  overridden_values_ = 0;
+  load_library_values();
+}
+
+AtomInfo::AtomInfo(const Ref<KeyVal>& keyval)
+{
+  initialize_names();
+  overridden_values_ = 0;
+  load_library_values();
+  override_library_values(keyval);
+}
+
+AtomInfo::AtomInfo(StateIn& s):
+  SavableState(s)
+{
+  initialize_names();
+  overridden_values_ = 0;
+  load_library_values();
+  char *overrides;
+  s.getstring(overrides);
+  if (overrides) {
+      Ref<ParsedKeyVal> keyval = new ParsedKeyVal;
+      keyval->parse_string(overrides);
+      override_library_values(keyval.pointer());
+      delete[] overrides;
+    }
+  if (s.version(::class_desc<AtomInfo>()) < 2) {
+      atomic_radius_scale_ = 1.0;
+      vdw_radius_scale_ = 1.0;
+      bragg_radius_scale_ = 1.0;
+      maxprob_radius_scale_ = 1.0;
+    }
+  else {
+      s.get(atomic_radius_scale_);
+      s.get(vdw_radius_scale_);
+      s.get(bragg_radius_scale_);
+      s.get(maxprob_radius_scale_);
+    }
+}
+
+AtomInfo::~AtomInfo()
+{
+  delete[] overridden_values_;
+}
+
+void
+AtomInfo::save_data_state(StateOut& s)
+{
+  s.putstring(overridden_values_);
+  s.put(atomic_radius_scale_);
+  s.put(vdw_radius_scale_);
+  s.put(bragg_radius_scale_);
+  s.put(maxprob_radius_scale_);
+}
+
+void
+AtomInfo::initialize_names()
+{
+  for (int i=0; i<Nelement; i++) {
+      symbol_to_Z_[elements_[i].symbol] = elements_[i].Z;
+      name_to_Z_[elements_[i].name] = elements_[i].Z;
+    }
+
+  // Z == DefaultZ is reserved for default values
+  symbol_to_Z_["Def"] = DefaultZ;
+  name_to_Z_["default"] = DefaultZ;
+
+  // Z == 1000 is reserved for point charges
+  name_to_Z_["charge"] = 1000;
+  symbol_to_Z_["Q"] = 1000;
+
+  for (std::map<std::string,int>::iterator i = symbol_to_Z_.begin();
+       i != symbol_to_Z_.end(); i++) {
+      Z_to_symbols_[i->second] = i->first;
+    }
+
+  for (std::map<std::string,int>::iterator i = name_to_Z_.begin();
+       i != name_to_Z_.end(); i++) {
+      Z_to_names_[i->second] = i->first;
+    }
+}
+
+void
+AtomInfo::load_library_values()
+{
+  Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
+  sc::auto_vec<char> in_char_array;
+  if (grp->me() == 0) {
+      const char* libdir;
+      std::string filename;
+      if ((libdir = getenv("SCLIBDIR")) != 0) {
+          const char* atominfo = "/atominfo.kv";
+          const char *eq = ::strchr(libdir,'=');
+          if (eq) libdir = eq + 1;
+          filename = std::string(libdir) + atominfo;
+        }
+      else {
+          struct stat sb;
+          const char *ainfo = SCDATADIR "/atominfo.kv";
+#ifdef SRC_SCLIBDIR
+          if (stat(ainfo, &sb) != 0) {
+              ainfo = SRC_SCLIBDIR "/atominfo.kv";
+            }
+#endif
+          filename = ainfo;
+        }
+      ExEnv::out0() << indent << "Reading file " << filename << "." << endl;
+      ifstream is(filename.c_str());
+      ostringstream ostrs;
+      is >> ostrs.rdbuf();
+      int n = 1 + strlen(ostrs.str().c_str());
+      in_char_array.reset(strcpy(new char[n],ostrs.str().c_str()));
+      grp->bcast(n);
+      grp->bcast(in_char_array.get(), n);
+    }
+  else {
+      int n;
+      grp->bcast(n);
+      in_char_array.reset(new char[n]);
+      grp->bcast(in_char_array.get(), n);
+    }
+
+  Ref<ParsedKeyVal> keyval = new ParsedKeyVal();
+  keyval->parse_string(in_char_array.get());
+  Ref<KeyVal> pkeyval = new PrefixKeyVal(keyval.pointer(), "atominfo");
+  load_values(pkeyval,0);
+}
+
+void
+AtomInfo::override_library_values(const Ref<KeyVal> &keyval)
+{
+  load_values(keyval, 1);
+}
+
+void
+AtomInfo::load_values(const Ref<KeyVal>& keyval, int override)
+{
+  Ref<Units> amu = new Units("amu");
+  Ref<Units> bohr = new Units("bohr");
+  Ref<Units> hartree = new Units("hartree");
+
+  load_values(Z_to_mass_, 0, "mass", keyval, override, amu);
+  load_values(Z_to_atomic_radius_, &atomic_radius_scale_, "atomic_radius",
+              keyval, override, bohr);
+  load_values(Z_to_vdw_radius_, &vdw_radius_scale_, "vdw_radius",
+              keyval, override, bohr);
+  load_values(Z_to_bragg_radius_, &bragg_radius_scale_,
+              "bragg_radius", keyval, override, bohr);
+  load_values(Z_to_maxprob_radius_, &maxprob_radius_scale_,
+              "maxprob_radius", keyval, override, bohr);
+  load_values(Z_to_rgb_, "rgb", keyval, override);
+  load_values(Z_to_ip_, 0, "ip", keyval, override, hartree);
+}
+
+void
+AtomInfo::load_values(std::map<int,double>&values,
+                      double *scale, const char *keyword,
+                      const Ref<KeyVal> &keyval, int override,
+                      const Ref<Units> &units)
+{
+  Ref<KeyVal> pkeyval = new PrefixKeyVal(keyval,keyword);
+  Ref<Units> fileunits = new Units(pkeyval->pcharvalue("unit"), Units::Steal);
+  double f = 1.0;
+  if (fileunits.nonnull() && units.nonnull()) {
+      f = fileunits->to(units);
+    }
+  double def = 0.0;
+  if (!override) {
+      def = pkeyval->doublevalue("default");
+      values[DefaultZ] = def;
+    }
+  int have_overridden = 0;
+  for (int elem=0; elem<Nelement; elem++) {
+      int Z = elements_[elem].Z;
+      double val = f * pkeyval->doublevalue(elements_[elem].symbol);
+      if (pkeyval->error() != KeyVal::OK) {
+          if (!override) values[Z] = def;
+        }
+      else {
+          values[Z] = val;
+          if (override) {
+              const char *prefix = " ";
+              if (!have_overridden) {
+                  add_overridden_value(keyword);
+                  add_overridden_value(":(");
+                  if (fileunits.nonnull() && fileunits->string_rep()) {
+                      char ustring[256];
+                      sprintf(ustring,"unit=\"%s\"",fileunits->string_rep());
+                      add_overridden_value(ustring);
+                    }
+                  else {
+                      prefix = "";
+                    }
+                  have_overridden = 1;
+                }
+              char *strval = pkeyval->pcharvalue(elements_[elem].symbol);
+              char assignment[256];
+              sprintf(assignment,"%s%s=%s", prefix,
+                      elements_[elem].symbol, strval);
+              delete[] strval;
+              add_overridden_value(assignment);
+            }
+        }
+    }
+  if (scale) {
+      KeyValValuedouble kvvscale(1.0);
+      *scale = pkeyval->doublevalue("scale_factor", kvvscale);
+      if (pkeyval->error() == KeyVal::OK) {
+          if (override) {
+              const char *prefix = " ";
+              if (!have_overridden) {
+                  add_overridden_value(keyword);
+                  add_overridden_value(":(");
+                  have_overridden = 1;
+                  prefix = "";
+                }
+              char *strval = pkeyval->pcharvalue("scale_factor");
+              char assignment[256];
+              sprintf(assignment,"%sscale_factor=%s", prefix, strval);
+              delete[] strval;
+              add_overridden_value(assignment);
+            }
+        }
+    }
+  if (have_overridden) {
+      add_overridden_value(")");
+    }
+}
+
+void
+AtomInfo::load_values(std::map<int,std::vector<double> >&values,
+                      const char *keyword,
+                      const Ref<KeyVal> &keyval, int override)
+{
+  Ref<KeyVal> pkeyval = new PrefixKeyVal(keyval,keyword);
+  double def[3];
+  if (!override) {
+      values[DefaultZ].resize(3);
+      for (int i=0; i<3; i++) {
+          def[i] = pkeyval->doublevalue("default",i);
+          values[DefaultZ][i] = def[i];
+        }
+    }
+  int have_overridden = 0;
+  for (int elem=0; elem<Nelement; elem++) {
+      double val;
+      int Z = elements_[elem].Z;
+      values[Z].resize(3);
+      for (int j=0; j<3; j++) {
+          val = pkeyval->doublevalue(elements_[elem].symbol,j);
+          if (pkeyval->error() != KeyVal::OK) {
+              if (!override) values[Z][j] = def[j];
+            }
+          else {
+              values[Z][j] = val;
+              if (override) {
+                  const char *prefix = " ";
+                  if (!have_overridden) {
+                      add_overridden_value(keyword);
+                      add_overridden_value(":(");
+                      prefix = "";
+                      have_overridden = 1;
+                    }
+                  char *strval = pkeyval->pcharvalue(elements_[elem].symbol,j);
+                  char assignment[256];
+                  sprintf(assignment,"%s%s:%d=%s",
+                          prefix, elements_[elem].symbol, j, strval);
+                  delete[] strval;
+                  add_overridden_value(assignment);
+                }
+            }
+        }
+    }
+  if (have_overridden) {
+      add_overridden_value(")");
+    }
+}
+
+void
+AtomInfo::add_overridden_value(const char *assignment)
+{
+  int length = strlen(assignment)+1;
+  if (overridden_values_) length += strlen(overridden_values_);
+  char *new_overridden_values = new char[length];
+  new_overridden_values[0] = '\0';
+  if (overridden_values_) strcat(new_overridden_values, overridden_values_);
+  strcat(new_overridden_values, assignment);
+  delete[] overridden_values_;
+  overridden_values_ = new_overridden_values;
+}
+
+int
+AtomInfo::string_to_Z(const std::string &name, int allow_exceptions)
+{
+  int Z;
+
+  // see if the name is a atomic number
+  Z = atoi(name.c_str());
+  if (Z) return Z;
+
+  // convert the name to lower case
+  std::string tmpname(name);
+  for (int j=0; j<tmpname.size(); j++) {
+      if (isupper(tmpname[j])) tmpname[j] = tolower(tmpname[j]);
+    }
+
+  std::map<std::string,int>::const_iterator iname
+      = name_to_Z_.find(tmpname);
+  if (iname != name_to_Z_.end()) return iname->second;
+
+  if (tmpname.size() > 0) {
+      if (islower(tmpname[0])) tmpname[0] = toupper(tmpname[0]);
+    }
+
+  iname = symbol_to_Z_.find(tmpname);
+  if (iname != symbol_to_Z_.end()) return iname->second;
+
+  if (allow_exceptions) {
+      ExEnv::err0() << sprintf("AtomInfo: invalid name: %s\n",name.c_str());
+      throw std::runtime_error("invalid atom name");
+    }
+
+  return 0;
+}
+
+double
+AtomInfo::lookup_value(const std::map<int,double>& values, int Z) const
+{
+  std::map<int,double>::const_iterator found = values.find(Z);
+  if (found == values.end()) {
+      found = values.find(DefaultZ);
+    }
+  return found->second;
+}
+
+double
+AtomInfo::lookup_array_value(const std::map<int,std::vector<double> >& values,
+                             int Z, int i) const
+{
+  std::map<int,std::vector<double> >::const_iterator found = values.find(Z);
+  if (found == values.end()) {
+      found = values.find(DefaultZ);
+    }
+  return found->second[i];
+}
+
+double
+AtomInfo::vdw_radius(int Z) const
+{
+  return lookup_value(Z_to_vdw_radius_,Z)*vdw_radius_scale_;
+}
+
+double
+AtomInfo::bragg_radius(int Z) const
+{
+  return lookup_value(Z_to_bragg_radius_,Z)*bragg_radius_scale_;
+}
+
+double
+AtomInfo::atomic_radius(int Z) const
+{
+  return lookup_value(Z_to_atomic_radius_,Z)*atomic_radius_scale_;
+}
+
+double
+AtomInfo::maxprob_radius(int Z) const
+{
+  return lookup_value(Z_to_maxprob_radius_,Z)*maxprob_radius_scale_;
+}
+
+double
+AtomInfo::ip(int Z) const
+{
+  return lookup_value(Z_to_ip_,Z);
+}
+
+double
+AtomInfo::rgb(int Z, int color) const
+{
+  return lookup_array_value(Z_to_rgb_,Z,color);
+}
+
+double
+AtomInfo::red(int Z) const
+{
+  return lookup_array_value(Z_to_rgb_,Z,0);
+}
+
+double
+AtomInfo::green(int Z) const
+{
+  return lookup_array_value(Z_to_rgb_,Z,1);
+}
+
+double
+AtomInfo::blue(int Z) const
+{
+  return lookup_array_value(Z_to_rgb_,Z,2);
+}
+
+double
+AtomInfo::mass(int Z) const
+{
+  return lookup_value(Z_to_mass_,Z);
+}
+
+std::string
+AtomInfo::name(int Z)
+{
+  std::map<int,std::string>::const_iterator found = Z_to_names_.find(Z);
+  if (found == Z_to_names_.end()) {
+      ExEnv::err0() << scprintf("AtomInfo: invalid Z: %d\n",Z);
+      throw std::runtime_error("invalid Z");
+    }
+  return found->second;
+}
+
+std::string
+AtomInfo::symbol(int Z)
+{
+  std::map<int,std::string>::const_iterator found = Z_to_symbols_.find(Z);
+  if (found == Z_to_symbols_.end()) {
+      ExEnv::err0() << scprintf("AtomInfo: invalid Z: %d\n",Z);
+      throw std::runtime_error("invalid Z");
+    }
+  return found->second;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
