shape.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

//
// shape.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[e-mail]>
// 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 <stdio.h>
#include <math.h>

#include <util/misc/formio.h>
#include <util/keyval/keyval.h>
#include <math/isosurf/shape.h>

using namespace std;
using namespace sc;

static const double shape_infinity = 1.0e23;

// given a vector X find which of the points in the vector of
// vectors, A, is closest to it and return the distance
static double
closest_distance(SCVector3& X,SCVector3*A,int n,SCVector3*grad)
{
  SCVector3 T = X-A[0];
  double min = T.dot(T);
  int imin = 0;
  for (int i=1; i<n; i++) {
      T = X-A[i];
      double tmp = T.dot(T);
      if (tmp < min) {min = tmp; imin = i;}
    }
  if (grad) {
      T = X - A[imin];
      T.normalize();
      *grad = T;
    }
  return sqrt(min);
}

//////////////////////////////////////////////////////////////////////
// Shape

static ClassDesc Shape_cd(
  typeid(Shape),"Shape",1,"public Volume",
  0, 0, 0);

Shape::Shape():
  Volume()
{
}

Shape::Shape(const Ref<KeyVal>& keyval):
  Volume(keyval)
{
}

Shape::~Shape()
{
}

void
Shape::compute()
{
  SCVector3 r;
  get_x(r);
  if (gradient_needed()) {
      if (!gradient_implemented()) {
          ExEnv::errn() << "Shape::compute: gradient not implemented" << endl;
          abort();
        }
      SCVector3 v;
      set_value(distance_to_surface(r,&v));
      set_actual_value_accuracy(desired_value_accuracy());
      set_gradient(v);
      set_actual_gradient_accuracy(desired_gradient_accuracy());
    }
  else if (value_needed()) {
      set_value(distance_to_surface(r));
      set_actual_value_accuracy(desired_value_accuracy());
    }
  if (hessian_needed()) {
      ExEnv::errn() << "Shape::compute(): can't do hessian yet" << endl;
      abort();
    }
}

int
Shape::is_outside(const SCVector3&r) const
{
  if (distance_to_surface(r)>0.0) return 1;
  return 0;
}

// Shape overrides volume's interpolate so it always gets points on
// the outside of the shape are always returned.

// interpolate using the bisection algorithm
void
Shape::interpolate(const SCVector3& A,
                   const SCVector3& B,
                   double val,
                   SCVector3& result)
{
  if (val < 0.0) {
      failure("Shape::interpolate(): val is < 0.0");
    }

  set_x(A);
  double value0 = value() - val;

  set_x(B);
  double value1 = value() - val;

  if (value0*value1 > 0.0) {
      failure("Shape::interpolate(): values at endpoints don't bracket val");
    }
  else if (value0 == 0.0) {
      result = A;
      return;
    }
  else if (value1 == 0.0) {
      result = B;
      return;
    }

  SCVector3 BA = B - A;

  double length = BA.norm();
  int niter = (int) (log(length/interpolation_accuracy())/M_LN2);

  double f0 = 0.0;
  double f1 = 1.0;
  double fnext = 0.5;

  SCVector3 X = A + fnext*BA;
  set_x(X);
  double valuenext = value() - val;

  do {
      for (int i=0; i<niter; i++) {
          if (valuenext*value0 <= 0.0) {
              value1 = valuenext;
              f1 = fnext;
              fnext = (f0 + fnext)*0.5;
            }
          else {
              value0 = valuenext;
              f0 = fnext;
              fnext = (fnext + f1)*0.5;
            }
          X = A + fnext*BA;
          set_x(X);
          valuenext = value() - val;
        }
      niter = 1;
    } while (valuenext < 0.0);

  result = X;
}

int
Shape::value_implemented() const
{
  return 1;
}

//////////////////////////////////////////////////////////////////////
// SphereShape

static ClassDesc SphereShape_cd(
  typeid(SphereShape),"SphereShape",1,"public Shape",
  0, create<SphereShape>, 0);

SphereShape::SphereShape(const SCVector3&o,double r):
  _origin(o),
  _radius(r)
{
}

SphereShape::SphereShape(const SphereShape&s):
  _origin(s._origin),
  _radius(s._radius)
{
}

SphereShape::SphereShape(const Ref<KeyVal>& keyval):
  _origin(new PrefixKeyVal(keyval,"origin")),
  _radius(keyval->doublevalue("radius"))
{
}

SphereShape::~SphereShape()
{
}

double
SphereShape::distance_to_surface(const SCVector3&p,SCVector3*grad) const
{
  int i;
  double r2 = 0.0;
  for (i=0; i<3; i++) {
      double tmp = p[i] - _origin[i];
      r2 += tmp*tmp;
    }
  double r = sqrt(r2);
  double d = r - _radius;
  if (grad) {
      SCVector3 pv(p);
      SCVector3 o(_origin);
      SCVector3 unit = pv - o;
      unit.normalize();
      for (i=0; i<3; i++) grad->elem(i) = unit[i];
    }
  return d;
}

void SphereShape::print(ostream&o) const
{
  o << indent
    << scprintf("SphereShape: r = %8.4f o = (%8.4f %8.4f %8.4f)",
                radius(),origin()[0],origin()[1],origin()[2])
    << endl;
}

void
SphereShape::boundingbox(double valuemin, double valuemax,
                         SCVector3& p1,
                         SCVector3& p2)
{
  if (valuemax < 0.0) valuemax = 0.0;

  int i;
  for (i=0; i<3; i++) {
      p1[i] = _origin[i] - _radius - valuemax;
      p2[i] = _origin[i] + _radius + valuemax;
    }
}

int
SphereShape::gradient_implemented() const
{
  return 1;
}

////////////////////////////////////////////////////////////////////////
// UncappedTorusHoleShape

static ClassDesc UncappedTorusHoleShape_cd(
  typeid(UncappedTorusHoleShape),"UncappedTorusHoleShape",1,"public Shape",
  0, 0, 0);

UncappedTorusHoleShape::UncappedTorusHoleShape(double r,
                               const SphereShape& s1,
                               const SphereShape& s2):
_s1(s1),
_s2(s2),
_r(r)
{
}

UncappedTorusHoleShape*
UncappedTorusHoleShape::newUncappedTorusHoleShape(double r,
                                                  const SphereShape&s1,
                                                  const SphereShape&s2)
{
  // if the probe sphere fits between the two spheres, then there
  // is no need to include this shape
  SCVector3 A(s1.origin());
  SCVector3 B(s2.origin());
  SCVector3 BA = B - A;
  if (2.0*r <  BA.norm() - s1.radius() - s2.radius()) return 0;

  // check to see if the surface is reentrant
  double rrs1 = r+s1.radius();
  double rrs2 = r+s2.radius();
  SCVector3 R12 = ((SCVector3)s1.origin()) - ((SCVector3)s2.origin());
  double r12 = sqrt(R12.dot(R12));
  if (sqrt(rrs1*rrs1-r*r) + sqrt(rrs2*rrs2-r*r) < r12)
    return new ReentrantUncappedTorusHoleShape(r,s1,s2);

  // otherwise create an ordinary torus hole
  return new NonreentrantUncappedTorusHoleShape(r,s1,s2);
}

// Given a node, finds a sphere in the plane of n and the centers
// of _s1 and _s2 that touches the UncappedTorusHole.  There are two
// candidates, the one closest to n is chosen.
void
UncappedTorusHoleShape::in_plane_sphere(
    const SCVector3& n,
    SCVector3& P) const
{
  // the center of the sphere is given by the vector equation
  // P = A + a R(AB) + b U(perp),
  // where
  // A is the center of _s1
  // B is the center of _s2
  // R(AB) is the vector from A to B, R(AB) = B - A
  // U(perp) is a unit vect perp to R(AB) and lies in the plane of n, A, and B
  // The unknown scalars, a and b are given by solving the following
  // equations:
  // | P - A | = r(A) + _r, and
  // | P - B | = r(B) + _r,
  // which give
  // | a R(AB) + b U(perp) | = r(A) + _r, and
  // | (a-1) R(AB) + b U(perp) | = r(B) + _r.
  // These further simplify to
  // a^2 r(AB)^2 + b^2 = (r(A)+_r)^2, and
  // (a-1)^2 r(AB)^2 + b^2 = (r(B)+_r)^2.
  // Thus,
  // a = (((r(A)+_r)^2 - (r(B)+_r)^2 )/(2 r(AB)^2)) + 1/2
  // b^2 = (r(A)+r)^2 - a^2 r(AB)^2

  SCVector3 A = _s1.origin();
  SCVector3 B = _s2.origin();
  SCVector3 N = n;
  SCVector3 R_AB = B - A;
  SCVector3 R_AN = N - A;

  // vector projection of R_AN onto R_AB
  SCVector3 P_AN_AB = R_AB * (R_AN.dot(R_AB)/R_AB.dot(R_AB));
  // the perpendicular vector
  SCVector3 U_perp = R_AN - P_AN_AB;

  // if |U| is tiny, then any vector perp to AB will do
  double u = U_perp.dot(U_perp);
  if (u<1.0e-23) {
      SCVector3 vtry = R_AB;
      vtry[0] += 1.0;
      vtry = vtry - R_AB * (vtry.dot(R_AB)/R_AB.dot(R_AB));
      if (vtry.dot(vtry) < 1.0e-23) {
          vtry = R_AB;
          vtry[1] += 1.0;
          vtry = vtry - R_AB * (vtry.dot(R_AB)/R_AB.dot(R_AB));
        }
      U_perp = vtry;
    }

  U_perp.normalize();
  //ExEnv::outn() << "A: "; A.print();
  //ExEnv::outn() << "U_perp: "; U_perp.print();
  //ExEnv::outn() << "R_AB: "; R_AB.print();

  double r_AB = sqrt(R_AB.dot(R_AB));
  double r_A = _s1.radius();
  double r_B = _s2.radius();

  double r_Ar = r_A + _r;
  double r_Br = r_B + _r;
  double a = ((r_Ar*r_Ar - r_Br*r_Br)/(2.0*r_AB*r_AB)) + 0.5;
  double b = sqrt(r_Ar*r_Ar - a*a*r_AB*r_AB);

  //ExEnv::outn() << scprintf("r_Ar = %f, r_AB = %f\n",r_Ar,r_AB);
  //ExEnv::outn() << scprintf("a = %f, b = %f\n",a,b);

  P = A + a * R_AB + b * U_perp;
  //ExEnv::outn() << "a*R_AB: "; (a*R_AB).print();
  //ExEnv::outn() << "b*U_perp: "; (b*U_perp).print();
}

void
UncappedTorusHoleShape::print(ostream&o) const
{
  o << indent << "UncappedTorusHoleShape:" << endl;
  o << incindent;
  o << indent << "r = " << _r << endl;
  o << indent << "s1 = ";
  o << incindent << skipnextindent;
  _s1.print(o);