Index: ThirdParty/mpqc_open/src/lib/chemistry/solvent/LIBS.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/solvent/LIBS.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/solvent/LIBS.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,2 @@
+libSCsolvent.LIBSUF
+#include <chemistry/molecule/LIBS.h>
Index: ThirdParty/mpqc_open/src/lib/chemistry/solvent/Makefile
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/solvent/Makefile	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/solvent/Makefile	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,72 @@
+#
+# Makefile
+#
+# Copyright (C) 1996 Limit Point Systems, Inc.
+#
+# Author: Curtis Janssen <cljanss@ca.sandia.gov>
+# 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.
+#
+
+TOPDIR=../../../..
+ifndef SRCDIR
+  SRCDIR=$(shell pwd)
+endif
+
+include $(SRCDIR)/$(TOPDIR)/lib/GlobalMakefile
+DEFINES += -DSRCDIR=\"$(SRCDIR)\"
+
+LD = $(CXX)
+CXXSRC = bem.cc disprep.cc
+
+LIBOBJ = $(CXXSRC:%.cc=%.$(OBJSUF))
+
+INC = bem.h
+
+DEPENDINCLUDE = $(INC)
+
+BIN_OR_LIB = LIB
+TARGET_TO_MAKE = libSCsolvent
+
+bemtest.testrun: DO_TESTRUN=no
+
+TESTSRC = bemtest.cc
+TESTOBJ = $(TESTSRC:%.cc=%.$(OBJSUF))
+TESTFILES = bemtest.in
+
+TESTPROGS = bemtest
+
+DISTFILES = $(CXXSRC) $(INC) Makefile $(TESTSRC) $(TESTFILES)
+
+default:: $(DEPENDINCLUDE)
+
+interface:: $(DEPENDINCLUDE)
+
+LIBS = $(shell $(LISTLIBS) $(INCLUDE) $(SRCDIR)/LIBS.h)
+
+bemtest: $(TESTOBJ) $(LIBS)
+	$(LTLINK) $(LD) $(LDFLAGS) -o bemtest $^ $(SYSLIBS) $(LTLINKBINOPTS)
+
+include $(SRCDIR)/$(TOPDIR)/lib/GlobalRules
+
+$(LIBOBJ:.$(OBJSUF)=.d): $(DEPENDINCLUDE)
+ifneq ($(DODEPEND),no)
+include $(LIBOBJ:.$(OBJSUF)=.d) $(TESTOBJ:%.$(OBJSUF)=%.d)
+endif
Index: ThirdParty/mpqc_open/src/lib/chemistry/solvent/bem.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/solvent/bem.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/solvent/bem.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,478 @@
+//
+// bem.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 <stdio.h>
+#include <util/misc/math.h>
+#include <util/misc/formio.h>
+#include <util/misc/timer.h>
+#include <math/scmat/matrix.h>
+#include <math/scmat/vector3.h>
+#include <math/scmat/local.h>
+#include <chemistry/solvent/bem.h>
+
+using namespace std;
+using namespace sc;
+
+static ClassDesc BEMSolvent_cd(
+  typeid(BEMSolvent),"BEMSolvent",1,"public DescribedClass",
+  0, create<BEMSolvent>, 0);
+
+BEMSolvent::BEMSolvent(const Ref<KeyVal>& keyval)
+{
+  vertex_area_ = 0;
+
+  matrixkit_ = new LocalSCMatrixKit;
+  
+  debug_ = keyval->intvalue("debug");
+  
+  solute_ << keyval->describedclassvalue("solute");
+
+  solvent_ << keyval->describedclassvalue("solvent");
+  // Use the aug-cc-pVQZ MP2 optimum geometry for H2O as default
+  if (solvent_.null()) {
+      solvent_ = new Molecule;
+      solvent_->add_atom(8, 0.0000000000,  0.0000000000, -0.1265941233);
+      solvent_->add_atom(1, 0.0000000000,  1.4304840085,  0.9856159541);
+      solvent_->add_atom(1, 0.0000000000, -1.4304840085,  0.9856159541);
+    }
+
+  solvent_density_ = keyval->doublevalue("solvent_density");
+  // use as default the number density of water in au^-3, T=25 C, P=101325 Pa
+  if (keyval->error() != KeyVal::OK) solvent_density_ = 0.004938887;
+
+  surf_ << keyval->describedclassvalue("surface");
+
+  dielectric_constant_ = keyval->doublevalue("dielectric_constant");
+  if (keyval->error() != KeyVal::OK) dielectric_constant_ = 78.0;
+
+  grp_ = MessageGrp::get_default_messagegrp();
+}
+
+BEMSolvent::~BEMSolvent()
+{
+}
+
+double**
+BEMSolvent::alloc_array(int n, int m)
+{
+  double ** result = new double*[n];
+  result[0] = new double[n*m];
+  for (int i=1; i<n; i++) {
+      result[i] = &result[i-1][m];
+    }
+  return result;
+}
+
+void
+BEMSolvent::free_array(double** array)
+{
+  if (!array) return;
+  delete[] array[0];
+  delete[] array;
+}
+
+void
+BEMSolvent::charge_positions(double**pos)
+{
+  int i,j;
+  int n = ncharge();
+  for (i=0; i<n; i++) {
+      const SCVector3& p = surf_->vertex(i)->point();
+      for (j=0; j<3; j++) { 
+          pos[i][j] = p[j];
+        }
+    }
+}
+
+void
+BEMSolvent::normals(double**norms)
+{
+  int i,j;
+  int n = ncharge();
+  for (i=0; i<n; i++) {
+      const SCVector3& p = surf_->vertex(i)->normal();
+      for (j=0; j<3; j++) { 
+          norms[i][j] = p[j];
+        }
+    }
+}
+
+void
+BEMSolvent::init()
+{
+  surf_->clear();
+  surf_->init();
+  system_matrix_i_ = 0;
+
+  f_ = (1.0-dielectric_constant_)/(2.0*M_PI*(1.0+dielectric_constant_));
+
+  if (vertex_area_) delete[] vertex_area_;
+  vertex_area_ = new double[ncharge()];
+  for (int i=0; i<ncharge(); i++) vertex_area_[i] = 0.0;
+  TriangulatedSurfaceIntegrator triint(surf_.pointer());
+  for (triint = 0; triint.update(); triint++) {
+      int j0 = triint.vertex_number(0);
+      int j1 = triint.vertex_number(1);
+      int j2 = triint.vertex_number(2);
+      double r = triint.r();
+      double s = triint.s();
+      double dA = triint.w();
+      vertex_area_[j0] += dA * (1 - r - s);
+      vertex_area_[j1] += dA * r;
+      vertex_area_[j2] += dA * s;
+    }
+}
+
+void
+BEMSolvent::done(int clear_surface)
+{
+  if (clear_surface) surf_->clear();
+  system_matrix_i_ = 0;
+
+  if (vertex_area_) delete[] vertex_area_;
+  vertex_area_ = 0;
+}
+
+void
+BEMSolvent::charges_to_surface_charge_density(double *charges)
+{
+  for (int i=0; i<ncharge(); i++) charges[i] /= vertex_area_[i];
+}
+
+void
+BEMSolvent::surface_charge_density_to_charges(double *charges)
+{
+  for (int i=0; i<ncharge(); i++) charges[i] *= vertex_area_[i];
+}
+
+double
+BEMSolvent::polarization_charge(double *charges)
+{
+  double charge = 0.0;
+  int n = ncharge();
+  for (int i=0; i<n; i++) charge += charges[i];
+  return charge;
+}
+
+// the passed enclosed_charge is determined by the called and
+// might different from the enclosed charge computed by Gauss's
+// law, which is stored as computed_enclosed_charge_
+void
+BEMSolvent::normalize_charge(double enclosed_charge, double* charges)
+{
+  int i;
+  double expected_charge = enclosed_charge
+                         * (1.0/dielectric_constant_ - 1.0);
+  double charge = 0.0;
+  double charge_pos = 0.0;
+  double charge_neg = 0.0;
+  int n = ncharge();
+  for (i=0; i<n; i++) {
+      charge += charges[i];
+      if (charges[i] > 0.0) charge_pos += charges[i];
+      else charge_neg += charges[i];
+    }
+
+  double scale_pos = 1.0;
+  double scale_neg = 1.0;
+  if (charge_pos > 1.0e-4 && charge_neg < -1.0e-4) {
+      scale_pos += (expected_charge-charge)/(2.0*charge_pos);
+      scale_neg += (expected_charge-charge)/(2.0*charge_neg);
+    }
+  else if (charge_pos > 1.0e-4) {
+      scale_pos += (expected_charge-charge)/charge_pos;
+    }
+  else if (charge_neg < -1.0e-4) {
+      scale_neg += (expected_charge-charge)/charge_neg;
+    }
+
+  double new_charge = 0.0;
+  for (i=0; i<n; i++) {
+      if (charges[i] > 0.0) charges[i] *= scale_pos;
+      else charges[i] *= scale_neg;
+      new_charge += charges[i];
+    }
+
+  if (fabs(new_charge - expected_charge) > 1.0e-3) {
+      ExEnv::outn() << "BEMSolvent:normalize_charge: failed:" << endl
+           << "new_charge = " << new_charge << endl
+           << "expected_charge = " << expected_charge << endl;
+      abort();
+    }
+
+  if (debug_) {
+      ExEnv::out0() << indent
+           << "BEMSolvent:normalize_charge:"
+           << endl << indent
+           << scprintf("  integrated surface charge = %20.15f", charge)
+           << endl << indent
+           << scprintf("  expected surface charge = %20.15f", expected_charge)
+           << endl;
+    }
+}
+
+void
+BEMSolvent::init_system_matrix()
+{
+  int i, j;
+  int n = ncharge();
+
+  RefSCDimension d = new SCDimension(n);
+  RefSCMatrix system_matrix(d,d,matrixkit());
+  system_matrix.assign(0.0);
+
+  tim_enter("precomp");
+  // precompute some arrays
+  TriangulatedSurfaceIntegrator triint(surf_.pointer());
+  int n_integration_points = triint.n();
+  SCVector3 *surfpv = new SCVector3[n_integration_points];
+  double *rfdA = new double[n_integration_points];
+  double *sfdA = new double[n_integration_points];
+  double *rsfdA = new double[n_integration_points];
+  int *j0 = new int[n_integration_points];
+  int *j1 = new int[n_integration_points];
+  int *j2 = new int[n_integration_points];
+  for (triint=0, i=0; i<n_integration_points&&triint.update(); i++,triint++) {
+      surfpv[i] = triint.current()->point();
+      j0[i] = triint.vertex_number(0);
+      j1[i] = triint.vertex_number(1);
+      j2[i] = triint.vertex_number(2);
+      double r = triint.r();
+      double s = triint.s();
+      double rs = 1 - r - s;
+      double dA = triint.w();
+      double fdA = - f_ * dA;
+      rfdA[i] = r * fdA;
+      sfdA[i] = s * fdA;
+      rsfdA[i] = rs * fdA;
+    }
+  tim_exit("precomp");
+
+  tim_enter("sysmat");
+  double *sysmati = new double[n];
+  RefSCVector vsysmati(system_matrix->rowdim(),system_matrix->kit());
+  // loop thru all the vertices
+  for (i = 0; i<n; i++) {
+      memset(sysmati,0,sizeof(double)*n);
+      Ref<Vertex> v = surf_->vertex(i);
+      const SCVector3& pv = v->point();
+      const SCVector3& nv = v->normal();
+      // integrate over the surface
+      for (j = 0; j < n_integration_points; j++) {
+          SCVector3 diff(pv - surfpv[j]);
+          double normal_component = diff.dot(nv);
+          double diff2 = diff.dot(diff);
+          if (diff2 <= 1.0e-8) {
+              // The self term must not be included here.  This
+              // case shouldn't occur for the usual integrators
+              // so abort.
+              ExEnv::errn() << "BEMSolvent: integrator gave the self term" << endl;
+              abort();
+            }
+          double denom = diff2*sqrt(diff2);
+          double common_factor = normal_component/denom;
+          sysmati[j0[j]] += common_factor * rsfdA[j];
+          sysmati[j1[j]] += common_factor * rfdA[j];
+          sysmati[j2[j]] += common_factor * sfdA[j];
+        }
+      vsysmati->assign(sysmati);
+      system_matrix->assign_row(vsysmati,i);
+    }
+  tim_exit("sysmat");
+
+  delete[] surfpv;
+  delete[] rfdA;
+  delete[] sfdA;
+  delete[] rsfdA;
+  delete[] j0;
+  delete[] j1;
+  delete[] j2;
+  delete[] sysmati;
+
+  tim_enter("AV");
+  double A = 0.0;
+  double V = 0.0;
+  for (triint = 0; triint.update(); triint++) {
+      V += triint.weight()*triint.dA()[2]*triint.current()->point()[2];
+      A += triint.w();
+    }
+  area_ = A;
+  volume_ = V;
+  tim_exit("AV");
+
+  ExEnv::out0() << indent
+       << scprintf("Solvent Accessible Surface:") << endl
+       << indent
+       << scprintf("  Area = %15.10f ", A)
+       << scprintf("Volume = %15.10f ", V)
+       << scprintf("Nvertex = %3d", n) << endl;
+
+  // Add I to the system matrix.
+  system_matrix->shift_diagonal(1.0);
+
+  //system_matrix->print("System Matrix");
+
+  tim_enter("inv");
+  system_matrix->invert_this();
+  system_matrix_i_ = system_matrix;
+  tim_exit("inv");
+
+  //system_matrix_i_->print("System Matrix Inverse");
+}
+
+void
+BEMSolvent::compute_charges(double* efield_dot_normals, double* charges)
+{
+  if (system_matrix_i_.null()) {
+      tim_enter("sysmat");
+      init_system_matrix();
+      tim_exit("sysmat");
+    }
+
+  tim_enter("qenq");
+  double efield_dot_normal = 0.0;
+  int n = ncharge();
+  for (int i=0; i<n; i++)
+      efield_dot_normal += efield_dot_normals[i] * vertex_area_[i];
+  tim_exit("qenq");
+
+  computed_enclosed_charge_ = efield_dot_normal/(4.0*M_PI);
+
+  if (debug_) {
+      double computed_expected_charge = computed_enclosed_charge_
+                                      * (1.0/dielectric_constant_ - 1.0);
+
+      ExEnv::out0() << indent
+         << scprintf("BEMSolvent:compute_charges: encl q = %20.15f",
+                     computed_enclosed_charge_)
+         << endl << indent
+         << scprintf("BEMSolvent:compute_charges: exp surface q = %20.15f",
+                     computed_expected_charge) << endl;
+    }
+
+  tim_enter("scomp");
+  RefSCVector edotn(system_matrix_i_.coldim(),matrixkit());
+  edotn.assign(efield_dot_normals);
+  //edotn.print("E dot normals");
+  edotn.scale(f_);
+  RefSCVector chrg = system_matrix_i_ * edotn;
+  //chrg.print("Charges");
+  chrg.convert(charges);
+  tim_exit("scomp");
+
+  tim_enter("stoq");
+  surface_charge_density_to_charges(charges);
+  tim_exit("stoq");
+}
+
+double
+BEMSolvent::nuclear_charge_interaction_energy(double *nuclear_charge,
+                                              double** charge_positions,
+                                              double* charge)
+{
+  double energy = 0.0;
+  int natom = solute_->natom();
+  for (int i=0; i<natom; i++) {
+      for (int j=0; j<ncharge(); j++) {
+          double r2 = 0.0;
+          for (int k=0; k<3; k++) {
+              double r = charge_positions[j][k] - solute_->r(i,k);
+              r2 += r*r;
+            }
+          energy += nuclear_charge[i] * charge[j] / sqrt(r2);
+        }
+    }
+  return energy;
+}
+
+double
+BEMSolvent::nuclear_interaction_energy(double** charge_positions,
+                                       double* charge)
+{
+  double energy = 0.0;
+  int natom = solute_->natom();
+  for (int i=0; i<natom; i++) {
+      for (int j=0; j<ncharge(); j++) {
+          double r2 = 0.0;
+          for (int k=0; k<3; k++) {
+              double r = charge_positions[j][k] - solute_->r(i,k);
+              r2 += r*r;
+            }
+          energy += double(solute_->Z(i)) * charge[j] / sqrt(r2);
+        }
+    }
+  return energy;
+}
+
+double
+BEMSolvent::self_interaction_energy(double** charge_positions,
+                                    double* charge)
+{
+  int i,j;
+
+  charges_to_surface_charge_density(charge);
+
+  TriangulatedSurfaceIntegrator triint(surf_.pointer());
+  int n_integration_points = triint.n();
+  SCVector3 *points = new SCVector3[n_integration_points];
+  double *charges = new double[n_integration_points];
+
+  double energy = 0.0;
+  for (triint=0, i=0; i<n_integration_points&&triint.update(); i++,triint++) {
+      points[i] = triint.current()->point();
+      int v0 = triint.vertex_number(0);
+      int v1 = triint.vertex_number(1);
+      int v2 = triint.vertex_number(2);
+      double r = triint.r();
+      double s = triint.s();
+      double rs = 1.0 - r - s;
+      double dA = triint.w();
+      charges[i] = (charge[v0]*rs + charge[v1]*r + charge[v2]*s)*dA;
+      energy += 0.0; // is this good enough for the self term?
+    }
+  for (i=0; i<n_integration_points; i++) {
+      double chargesi = charges[i];
+      SCVector3 pointsi(points[i]);
+      for (j = 0; j<i; j++) {
+          energy += chargesi*charges[j]/pointsi.dist(points[j]);
+        }
+    }
+
+  delete[] points;
+  delete[] charges;
+
+  surface_charge_density_to_charges(charge);
+
+  return energy;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/chemistry/solvent/bem.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/solvent/bem.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/solvent/bem.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,161 @@
+//
+// bem.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_solvent_bem_h
+#define _chemistry_solvent_bem_h
+
+#include <util/class/class.h>
+#include <util/state/state.h>
+#include <util/keyval/keyval.h>
+#include <math/isosurf/volume.h>
+#include <math/isosurf/surf.h>
+#include <math/scmat/matrix.h>
+#include <chemistry/molecule/molecule.h>
+
+namespace sc {
+
+// This represents a solvent by a polarization charge on a dielectric
+// boundary surface.
+class BEMSolvent: public DescribedClass {
+  private:
+    int debug_;
+
+    Ref<Molecule> solute_;
+    Ref<Molecule> solvent_;
+    double solvent_density_;
+    double dielectric_constant_;
+    Ref<SCMatrixKit> matrixkit_;
+    RefSCMatrix system_matrix_i_;
+    double f_;
+    Ref<MessageGrp> grp_;
+
+    double area_;
+    double volume_;
+    double computed_enclosed_charge_;
+    double edisp_;
+    double erep_;
+
+    Ref<TriangulatedImplicitSurface> surf_;
+
+    double** alloc_array(int n, int m);
+    void free_array(double**);
+
+    // This holds the area associated with each vertex.  It is used
+    // to convert charges to charge densities and back.
+    double* vertex_area_;
+
+    // Given charges compute surface charge density.
+    void charges_to_surface_charge_density(double *charges);
+
+    // Given surface charge density compute charges.
+    void surface_charge_density_to_charges(double *charges);
+  public:
+    BEMSolvent(const Ref<KeyVal>&);
+    virtual ~BEMSolvent();
+
+    // This should be called after everything is setup--the
+    // molecule has the correct the geometry and all of the
+    // parameters have been adjusted.
+    void init();
+    // This gets rid of the system matrix inverse and, optionally, the surface.
+    void done(int clear_surface = 1);
+
+    int ncharge() { return surf_->nvertex(); }
+
+    Ref<Molecule> solvent() { return solvent_ ;}
+    double solvent_density() { return solvent_density_ ;}
+
+    // NOTE: call allocation routines after init and free routines before done
+    double** alloc_charge_positions() { return alloc_array(ncharge(), 3); }
+    void free_charge_positions(double**a) { free_array(a); }
+
+    double** alloc_normals()  { return alloc_array(ncharge(), 3); }
+    void free_normals(double**a) { free_array(a); }
+
+    double* alloc_efield_dot_normals()  { return new double[ncharge()]; }
+    void free_efield_dot_normals(double*a) { delete[] a; }
+
+    double* alloc_charges() { return new double[ncharge()]; }
+    void free_charges(double*a) { delete[] a; }
+
+    void charge_positions(double**);
+    void normals(double**);
+
+    // Given the efield dotted with the normals at the charge positions this
+    // will compute a new set of charges.
+    void compute_charges(double* efield_dot_normals, double* charge);
+
+    // Given a set of charges and a total charge, this will normalize
+    // the integrated charge to the charge that would be expected on
+    // the surface if the given total charge were enclosed within it.
+    void normalize_charge(double enclosed_charge, double* charges);
+
+    // Given charges and nuclear charges compute their interation energy.
+    double nuclear_charge_interaction_energy(double *nuclear_charge,
+                                             double** charge_positions,
+                                             double* charge);
+
+    // Given charges compute the interaction energy between the nuclei
+    // and the point charges.
+    double nuclear_interaction_energy(double** charge_positions,
+                                      double* charge);
+
+    // Given charges compute the interaction energy for just the surface.
+    double self_interaction_energy(double** charge_positions, double *charge);
+    
+    // Given the charges, return the total polarization charge on the surface.
+    double polarization_charge(double* charge);
+
+    // Return the area (available after compute_charges called).
+    double area() const { return area_; }
+    // Return the volume (available after compute_charges called).
+    double volume() const { return volume_; }
+    // Return the enclosed charge (available after compute_charges called).
+    double computed_enclosed_charge() const {
+      return computed_enclosed_charge_;
+    }
+
+    double disp() {return edisp_;}
+    double rep()  {return erep_;}
+    double disprep();
+
+    // this never needs to be called explicitly, but is here now for debugging
+    void init_system_matrix();
+
+    Ref<TriangulatedImplicitSurface> surface() const { return surf_; }
+
+    Ref<SCMatrixKit> matrixkit() { return matrixkit_; }
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/chemistry/solvent/bemtest.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/solvent/bemtest.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/solvent/bemtest.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,66 @@
+//
+// bemtest.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 <stdio.h>
+#include <util/misc/ieee.h>
+#include <util/keyval/keyval.h>
+#include <chemistry/molecule/molshape.h>
+#include <chemistry/solvent/bem.h>
+
+using namespace sc;
+
+int
+main()
+{
+  // Abort on floating point errors.
+  ieee_trap_errors();
+
+  Ref<KeyVal> keyval = new ParsedKeyVal(SRCDIR "/bemtest.in");
+
+  Ref<BEMSolvent> solvent;
+  solvent << keyval->describedclassvalue("solvent");
+
+  solvent->init();
+  solvent->init_system_matrix();
+  solvent->done();
+
+  solvent->init();
+  solvent->init_system_matrix();
+  solvent->done();
+
+  ConnollyShape::print_counts();
+  CS2Sphere::print_counts();
+
+  return 0;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/chemistry/solvent/bemtest.in
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/solvent/bemtest.in	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/solvent/bemtest.in	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,40 @@
+% Emacs should use -*- keyval -*- mode.
+
+h2o<Molecule>: (
+    symmetry=c2v
+    {atoms geometry} = {
+      H  [   1.5  0.0  -0.3 ]
+      O  [   0.0  0.0   1.0 ]
+     }
+  )
+
+he<Molecule>: (
+    symmetry=c1
+    {atoms geometry} = {
+      He  [   0.0  0.0  0.0 ]
+     }
+  )
+
+molecule = $:h2o
+
+atominfo<AtomInfo>: ()
+
+solvent<BEMSolvent>: (
+    solute = $:molecule
+    surface<TriangulatedImplicitSurface>: (
+        surface<TriangulatedSurface>: (
+            verbose = yes
+          )
+        volume<ConnollyShape>: (
+            molecule = $:molecule
+            probe_radius = 2.6456173
+            atominfo = $:atominfo
+          )
+        resolution = 0.5
+        remove_short_edges = 1
+        remove_slender_triangles = 1
+        short_edge_factor = 0.4
+        slender_triangle_factor = 0.2
+      )
+    dielectric_constant = 80.0
+  )
Index: ThirdParty/mpqc_open/src/lib/chemistry/solvent/disprep.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/chemistry/solvent/disprep.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/chemistry/solvent/disprep.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,335 @@
+
+#ifdef HAVE_CONFIG_H
+#include <scconfig.h>
+#endif
+
+#include <fstream>
+
+#include <util/keyval/keyval.h>
+#include <math/isosurf/shape.h>
+#include <chemistry/qc/wfn/solvent.h>
+#include <chemistry/molecule/formula.h>
+
+#ifdef USING_NAMESPACE_STD
+using namespace std;
+#endif
+using namespace sc;
+
+static inline double
+get_ki(int z)
+{
+  // The ki values (used in the computation of the dispersion coefficients)
+  // for H, C, and N were taken from Vigne-Maeder and Claverie, JACS 1987, v109, pp24-28
+  // and the value for O from Huron and Claverie, J. Phys. Chem. 1974, v78, p1862
+
+  double ki;
+  
+  if (z <= 0) {
+      ExEnv::errn() << "Non-positive nuclear charge encountered in computation of"
+           << " dispersion coefficient" << endl;
+      abort();
+    }
+  else if (z == 1) ki = 1.0;
+  else if (z == 6) ki = 1.0;
+  else if (z == 7) ki = 1.18;
+  else if (z == 8) ki = 1.36;  // from Huron & Claverie, J.Phys.Chem v78, 1974, p1862
+  else if (z > 1 && z < 6) {
+      ki = 1.0;
+      ExEnv::out0() << "Warning: No d6 dispersion coefficient available for atomic number " <<
+              z << "; using value for carbon instead" << endl;
+    }
+  else {
+      ki = 1.18;
+      ExEnv::out0() << "Warning: No d6 dispersion coefficient available for atomic number " <<
+              z << "; using value for nitrogen instead" << endl;
+    }
+  
+  return ki;
+}
+
+static inline double
+get_d6ii(int z, double r_vdw)
+{
+  // The dispersion coefficient d6 for a pair of atoms ij can be computed
+  // from the dispersion coefficient d6ii for atom pair ii and d6jj for
+  // atom pair jj by the formula: d6 = sqrt(d6ii*d6jj).
+  // The dispersion coefficients d8 and d10 can be obtained from d6.
+  // The d6ii values given below were taken from: Vigne-Maeder and Claverie
+  // JACS 1987, v. 109, pp. 24-28.
+
+  const double a6 = 0.143; // [kcal/mol]
+  double d6ii;
+  double ki;
+
+  Ref<Units> unit = new Units("kcal/mol");
+  
+  ki = get_ki(z);
+  d6ii = ki*ki*a6*pow(4*r_vdw*r_vdw,3.0);  // units of (kcal mol^-1)*bohr^6
+  d6ii *= unit->to_atomic_units();  // convert to atomic units
+  return d6ii;
+}
+
+static inline double
+get_d8ii(double d6ii, double r_vdw)
+{
+  // The value of c8 was taken from Vigne-Maeder and Claverie, JACS 1987,
+  // v. 109, pp 24-28 and is here obtained in atomic units by using
+  // atomic units for d6ii and r_vdw
+
+  double d8ii;
+  const double c8 = 0.26626;
+
+  d8ii = d6ii*c8*4*pow(r_vdw,2.0);
+  
+  return d8ii;
+}
+
+static inline double
+get_d10ii(double d6ii, double r_vdw)
+{
+  // The value of c10 was taken from Vigne-Maeder and Claverie, JACS 1987,
+  // v. 109, pp 24-28 and is here obtained in atomic units by using
+  // atomic units for d6ii and r_vdw
+
+  double d10ii;
+  const double c10 = 0.095467;
+
+  d10ii = d6ii*c10*16*pow(r_vdw,4.0);
+  
+  return d10ii;
+}
+
+// For debugging compute 6, 8, and 10 contributions separately
+static inline double
+disp6_contrib(double rasnorm, double d6)
+{
+  double edisp6_contrib;
+
+  edisp6_contrib = d6/(3*pow(rasnorm,6.0)); // atomic units
+  
+  return edisp6_contrib;
+}
+
+static inline double
+disp8_contrib(double rasnorm, double d8)
+{
+  double edisp8_contrib;
+
+  edisp8_contrib = d8/(5*pow(rasnorm,8.0)); // atomic units
+  
+  return edisp8_contrib;
+}
+
+static inline double
+disp10_contrib(double rasnorm, double d10)
+{
+  double edisp10_contrib;
+
+  edisp10_contrib = d10/(7*pow(rasnorm,10.0)); // atomic units
+  
+  return edisp10_contrib;
+}
+
+static inline double
+disp_contrib(double rasnorm, double d6, double d8, double d10)
+{
+  double edisp_contrib;
+
+  edisp_contrib = d6/(3*pow(rasnorm,6.0)) + d8/(5*pow(rasnorm,8.0))
+                + d10/(7*pow(rasnorm,10.0));
+  
+  return edisp_contrib;
+}
+
+static inline double
+rep_contrib(double rasnorm, double ri_vdw, double rj_vdw, double ki, double kj,
+            double kcalpermol_to_hartree)
+{
+  // The expression and the parameters used for the repulsion energy
+  // were taken from Vigne-Maeder and Claverie, JACS 1987, v109, pp24-28
+  // NB: We have omitted the factor Gij
+  
+  const double c = 90000; // [kcal/mol]
+  const double gamma = 12.35;
+  double erep_contrib;
+  double tmp;
+  
+  tmp = gamma*rasnorm/(2.0*sqrt(ri_vdw*rj_vdw));
+
+  erep_contrib = -ki*kj*c*(1.0/tmp + 2.0/(tmp*tmp) + 2.0/(tmp*tmp*tmp))*exp(-tmp);
+  erep_contrib *= kcalpermol_to_hartree; // convert from kcal/mol to atomic units
+  
+  return erep_contrib;
+}
+
+double
+BEMSolvent::disprep() 
+{
+  double edisprep = 0.0;
+  double edisprep_contrib;
+  double edisp6_contrib, edisp8_contrib, edisp10_contrib; // for debugging
+  double erep_contrib;
+  double edisp6 = 0.0; // for debugging
+  double edisp8 = 0.0; // for debugging
+  double edisp10 = 0.0; // for debugging
+  double erep = 0.0;
+  double proberadius;
+  double radius;
+  double rasnorm;
+  double weight;
+  double d6, d8, d10; // dispersion coefficients
+  double d6aa, d8aa, d10aa; // dispersion coefficients for atom pair aa
+  double d6ss, d8ss, d10ss; // dispersion coefficients for atom pair ss
+  int i, iloop, isolute;
+  int natomtypes;
+  int z_solvent_atom;
+
+  Ref<Units> unit = new Units("kcal/mol");
+  double kcalpermol_to_hartree = unit->to_atomic_units();
+
+  Ref<AtomInfo> atominfo = solute_->atominfo();
+  Ref<AtomInfo> solventatominfo = solvent_->atominfo();
+  MolecularFormula formula(solvent_);
+
+  // Compute number of different atom types in solvent molecule
+  natomtypes = formula.natomtypes();
+
+  double *solute_d6ii  = new double[solute_->natom()];
+  double *solute_d8ii  = new double[solute_->natom()];
+  double *solute_d10ii = new double[solute_->natom()];
+  double *solute_ki = new double[solute_->natom()];
+
+  for (isolute=0; isolute<solute_->natom(); isolute++) {
+      int Z_solute = solute_->Z(isolute);
+      double radius = atominfo->vdw_radius(Z_solute);
+      solute_d6ii[isolute] = get_d6ii(Z_solute,radius);
+      solute_d8ii[isolute] = get_d8ii(solute_d6ii[isolute],radius);
+      solute_d10ii[isolute] = get_d10ii(solute_d6ii[isolute],radius);
+      solute_ki[isolute] = get_ki(Z_solute);
+    }
+  
+  // Loop over atom types in solvent molecule
+  for (iloop=0; iloop<natomtypes; iloop++) {
+
+      // define the shape of the surface for current atom type
+      Ref<UnionShape> us = new UnionShape;
+      z_solvent_atom = formula.Z(iloop);
+      proberadius = solventatominfo->vdw_radius(z_solvent_atom);
+      for (i=0; i<solute_->natom(); i++) {
+          us->add_shape(new SphereShape(solute_->r(i),
+                                        atominfo->vdw_radius(solute_->Z(i))+proberadius));
+        }
+      
+      // triangulate the surface
+      Ref<AssignedKeyVal> keyval = new AssignedKeyVal;
+      keyval->assign("volume", us.pointer());
+      keyval->assign("order", 2);
+      keyval->assign("remove_short_edges", 1);
+      keyval->assign("remove_small_triangles", 1);
+      keyval->assign("remove_slender_triangles", 1);
+      keyval->assign("short_edge_factor", 0.8);
+      keyval->assign("small_triangle_factor", 0.8);
+      keyval->assign("slender_triangle_factor", 0.8);
+      Ref<TriangulatedImplicitSurface> ts = new TriangulatedImplicitSurface(keyval.pointer());
+      ts->init();
+
+      // Debug print: check the triangulated surface
+//      if (iloop == 0) {
+//          ofstream geomviewfile("geomview.input");
+//          ts->print_geomview_format(geomviewfile);
+//        }
+      
+      ExEnv::out0().setf(ios::scientific,ios::floatfield); // use scientific format
+      ExEnv::out0() << "Area of disp-rep surface generated with atom number "
+           << setw(3) << setfill(' ') << z_solvent_atom
+           << " as probe: " << setprecision(4) << ts->area()
+           << " bohr^2" << endl;
+      
+      edisprep_contrib = 0.0;
+      edisp6_contrib = 0.0;  // for debugging
+      edisp8_contrib = 0.0;  // for debugging
+      edisp10_contrib = 0.0; // for debugging
+      erep_contrib = 0.0;
+      TriangulatedSurfaceIntegrator triint(ts.pointer());
+
+      double solvent_ki = get_ki(z_solvent_atom);
+      d6ss = get_d6ii(z_solvent_atom,proberadius);
+      d8ss = get_d8ii(d6ss, proberadius);
+      d10ss = get_d10ii(d6ss, proberadius);
+              
+      // integrate the surface
+      for (triint=0; triint.update(); triint++) {
+          SCVector3 dA = triint.dA();
+          SCVector3 location = triint.current()->point();
+          weight = triint.weight();
+          
+          //Loop over atoms in solute
+          for (isolute=0; isolute<solute_->natom(); isolute++) {
+
+              SCVector3 atom(solute_->r(isolute)); 
+              SCVector3 ras = location - atom;
+              rasnorm = ras.norm();
+              radius = atominfo->vdw_radius(solute_->Z(isolute));
+              d6aa = solute_d6ii[isolute];
+              d8aa = solute_d8ii[isolute];
+              d10aa = solute_d10ii[isolute];
+              d6 = sqrt(d6aa*d6ss);
+              d8 = sqrt(d8aa*d8ss);
+              d10 = sqrt(d10aa*d10ss);
+
+              double f = ras.dot(dA)*weight;
+              double tdisp6 = f*disp6_contrib(rasnorm,d6);
+              double tdisp8 = f*disp8_contrib(rasnorm,d8);
+              double tdisp10 = f*disp10_contrib(rasnorm,d10);
+              double trep = f*rep_contrib(rasnorm,radius,proberadius,
+                                          solute_ki[isolute],solvent_ki,
+                                          kcalpermol_to_hartree);
+              double tdisp = tdisp6+tdisp8+tdisp10;
+
+              // add in contributions to various energies; the minus sign
+              // is there to get the normal pointing into the cavity
+              edisprep_contrib -= tdisp+trep;
+              edisp6_contrib -= tdisp6;
+              edisp8_contrib -= tdisp8;
+              edisp10_contrib -= tdisp10;
+              erep_contrib -= trep;
+              
+            }
+        }
+      
+      edisprep += edisprep_contrib*formula.nZ(iloop);
+      edisp6 += edisp6_contrib*formula.nZ(iloop);
+      edisp8 += edisp8_contrib*formula.nZ(iloop);
+      edisp10 += edisp10_contrib*formula.nZ(iloop);
+      erep += erep_contrib*formula.nZ(iloop);
+    }
+
+  delete[] solute_d6ii;
+  delete[] solute_d8ii;
+  delete[] solute_d10ii;
+  delete[] solute_ki;
+
+  // Multiply energies by number density of solvent
+  // Print out individual energy contributions in kcal/mol
+  
+  ExEnv::out0().setf(ios::scientific,ios::floatfield); // use scientific format
+  ExEnv::out0().precision(5);
+  ExEnv::out0() << "Edisp6:  " << edisp6*solvent_density_*unit->from_atomic_units()
+       << " kcal/mol" << endl; 
+  ExEnv::out0() << "Edisp8:  " << edisp8*solvent_density_*unit->from_atomic_units()
+       << " kcal/mol" << endl;
+  ExEnv::out0() << "Edisp10: " << edisp10*solvent_density_*unit->from_atomic_units()
+       << " kcal/mol" << endl;
+
+
+  ExEnv::out0() << "Total dispersion energy: "
+       << (edisp6 + edisp8 + edisp10)*solvent_density_*unit->from_atomic_units()
+       << " kcal/mol" << endl;
+  ExEnv::out0() << "Repulsion energy:        " << setw(12) << setfill(' ') 
+       << erep*solvent_density_*unit->from_atomic_units() << " kcal/mol" << endl;
+  
+  return edisprep*solvent_density_; // atomic units
+
+}
+
+
