shape.cc

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

                  }
              }
          }
        else _intersects_AB = 0;
      }

      {
        // Does the circle of intersection intersect with BC?
        SCVector3 MC = M - C();
        SCVector3 uBC = B() - C(); uBC.normalize();
        SCVector3 XC = uBC * MC.dot(uBC);
        SCVector3 XM = XC - MC;
        double rXM2 = XM.dot(XM);
        double d_intersect_from_x2 = r_intersect*r_intersect - rXM2;
        if (d_intersect_from_x2 > 0.0) {
            _intersects_BC = 1;
            double tmp = sqrt(d_intersect_from_x2);
            double d_intersect_from_x[2];
            d_intersect_from_x[0] = tmp;
            d_intersect_from_x[1] = -tmp;
            for (i=0; i<2; i++) {
                for (int j=0; j<2; j++) {
                    IBCD[i][j] = XM + d_intersect_from_x[j]*uBC + MD[i];
                  }
              }
          }
        else _intersects_BC = 0;
      }

      {
        // Does the circle of intersection intersect with CA?
        SCVector3 MA = M - A();
        SCVector3 uCA = C() - A(); uCA.normalize();
        SCVector3 XA = uCA * MA.dot(uCA);
        SCVector3 XM = XA - MA;
        double rXM2 = XM.dot(XM);
        double d_intersect_from_x2 = r_intersect*r_intersect - rXM2;
        if (d_intersect_from_x2 > 0.0) {
            _intersects_CA = 1;
            double tmp = sqrt(d_intersect_from_x2);
            double d_intersect_from_x[2];
            d_intersect_from_x[0] = tmp;
            d_intersect_from_x[1] = -tmp;
            for (i=0; i<2; i++) {
                for (int j=0; j<2; j++) {
                    ICAD[i][j] = XM + d_intersect_from_x[j]*uCA + MD[i];
                  }
              }
          }
        else _intersects_CA = 0;
      }

    }

#if 0 // test code
  ExEnv::outn() << "Uncapped5SphereExclusionShape: running some tests" << endl;
  verbose = 1;

  FILE* testout = fopen("testout.vect", "w");

  const double scalefactor_inc = 0.05;
  const double start = -10.0;
  const double end = 10.0;

  SCVector3 middle = (1.0/3.0)*(A()+B()+C());

  int nlines = 1;
  int nvert = (int) ( (end-start) / scalefactor_inc);
  int ncolor = nvert;

  fprintf(testout, "VECT\n%d %d %d\n", nlines, nvert, ncolor);

  fprintf(testout, "%d\n", nvert);
  fprintf(testout, "%d\n", ncolor);

  double scalefactor = start;
  for (int ii = 0; ii<nvert; ii++) {
      SCVector3 position = (D[0] - middle) * scalefactor + middle;
      double d = distance_to_surface(position);
      fprintf(testout, "%f %f %f # value = %f\n",
              position[0], position[1], position[2], d);
      scalefactor += scalefactor_inc;
    }
  scalefactor = start;
  for (ii = 0; ii<nvert; ii++) {
      SCVector3 position = (D[0] - middle) * scalefactor + middle;
      double d = distance_to_surface(position);
      ExEnv::outn() << scprintf("d = %f\n", d);
      if (d<0.0) fprintf(testout,"1.0 0.0 0.0 0.5\n");
      else fprintf(testout,"0.0 0.0 1.0 0.5\n");
      scalefactor += scalefactor_inc;
    }

  fclose(testout);
  ExEnv::outn() << "testout.vect written" << endl;

  verbose = 0;
#endif // test code
  
}
int
Uncapped5SphereExclusionShape::is_outside(const SCVector3&X) const
{
  SCVector3 Xv(X);

  if (verbose) ExEnv::outn() << scprintf("point %14.8f %14.8f %14.8f\n",X(0),X(1),X(2));

  // The folded case isn't handled correctly here, so use
  // the less efficient distance_to_surface routine.
  if (_folded) {
      return distance_to_surface(X) >= 0.0;
    }

  for (int i=0; i<2; i++) {
      SCVector3 XD = Xv - D[i];
      double rXD = XD.norm();
      if (rXD <= r()) return 1;
      double u = BDxCD[i].dot(XD)/BDxCDdotAD[i];
      if (u <= 0.0) return 1;
      double v = CDxAD[i].dot(XD)/CDxADdotBD[i];
      if (v <= 0.0) return 1;
      double w = ADxBD[i].dot(XD)/ADxBDdotCD[i];
      if (w <= 0.0) return 1;
    }

  if (verbose) ExEnv::outn() << "is_inside" << endl;

  return 0;
}
static int
is_contained_in_unbounded_pyramid(SCVector3 XD,
                                  SCVector3 AD,
                                  SCVector3 BD,
                                  SCVector3 CD)
{
  SCVector3 BDxCD = BD.cross(CD);
  SCVector3 CDxAD = CD.cross(AD);
  SCVector3 ADxBD = AD.cross(BD);
  double u = BDxCD.dot(XD)/BDxCD.dot(AD);
  if (u <= 0.0) return 0;
  double v = CDxAD.dot(XD)/CDxAD.dot(BD);
  if (v <= 0.0) return 0;
  double w = ADxBD.dot(XD)/ADxBD.dot(CD);
  if (w <= 0.0) return 0;
  return 1;
}
double
Uncapped5SphereExclusionShape::
  distance_to_surface(const SCVector3&X,SCVector3*grad) const
{
  SCVector3 Xv(X);

  // Find out if I'm on the D[0] side or the D[1] side of the ABC plane.
  int side;
  SCVector3 XM = Xv - M;
  if (MD[0].dot(XM) > 0.0) side = 1;
  else side = 0;

  if (verbose) {
      ExEnv::outn() << scprintf("distance_to_surface: folded = %d, side = %d\n",
                       _folded, side);
      ExEnv::outn() << "XM = "; XM.print();
      ExEnv::outn() << "MD[0] = "; MD[0].print();
      ExEnv::outn() << "MD[0].dot(XM) = " << MD[0].dot(XM) << endl;
    }

  SCVector3 XD = Xv - D[side];
  double u = BDxCD[side].dot(XD)/BDxCDdotAD[side];
  if (verbose) ExEnv::outn() << scprintf("u = %14.8f\n", u);
  if (u <= 0.0) return shape_infinity;
  double v = CDxAD[side].dot(XD)/CDxADdotBD[side];
  if (verbose) ExEnv::outn() << scprintf("v = %14.8f\n", v);
  if (v <= 0.0) return shape_infinity;
  double w = ADxBD[side].dot(XD)/ADxBDdotCD[side];
  if (verbose) ExEnv::outn() << scprintf("w = %14.8f\n", w);
  if (w <= 0.0) return shape_infinity;
  double rXD = XD.norm();
  if (verbose) ExEnv::outn() << scprintf("r() - rXD = %14.8f\n", r() - rXD);
  if (rXD <= r()) {
      if (!_reentrant) return r() - rXD;
      // this shape is reentrant
      // make sure that we're on the right side
      if ((side == 1) || (u + v + w <= 1.0)) {
          // see if we're outside the cone that intersects
          // the opposite sphere
          double angle = acos(fabs(XD.dot(MD[side]))
                              /(XD.norm()*MD[side].norm()));
          if (angle >= theta_intersect) {
              if (grad) {
                  *grad = (-1.0/rXD)*XD;
                }
              return r() - rXD;
            }
          if (_intersects_AB
              &&is_contained_in_unbounded_pyramid(XD,
                                                  MD[side],
                                                  IABD[side][0],
                                                  IABD[side][1])) {
              //ExEnv::outn() << scprintf("XD: "); XD.print();
              //ExEnv::outn() << scprintf("MD[%d]: ",i); MD[i].print();
              //ExEnv::outn() << scprintf("IABD[%d][0]: ",i); IABD[i][0].print();
              //ExEnv::outn() << scprintf("IABD[%d][1]: ",i); IABD[i][1].print();
              return closest_distance(XD,(SCVector3*)IABD[side],2,grad);
            }
          if (_intersects_BC
              &&is_contained_in_unbounded_pyramid(XD,
                                                  MD[side],
                                                  IBCD[side][0],
                                                  IBCD[side][1])) {
              return closest_distance(XD,(SCVector3*)IBCD[side],2,grad);
            }
          if (_intersects_CA
              &&is_contained_in_unbounded_pyramid(XD,
                                                  MD[side],
                                                  ICAD[side][0],
                                                  ICAD[side][1])) {
              return closest_distance(XD,(SCVector3*)ICAD[side],2,grad);
            }
          // at this point we are closest to the ring formed
          // by the intersection of the two probe spheres
          double distance_to_plane;
          double distance_to_ring_in_plane;
          double MDnorm = MD[side].norm();
          SCVector3 XM = XD - MD[side];
          SCVector3 XM_in_plane;
          if (MDnorm<1.0e-16) {
              distance_to_plane = 0.0;
              XM_in_plane = XD;
            }
          else {
              distance_to_plane = XM.dot(MD[side])/MDnorm;
              XM_in_plane = XM - (distance_to_plane/MDnorm)*MD[side];
            }
          if (grad) {
              double XM_in_plane_norm = XM_in_plane.norm();
              if (XM_in_plane_norm < 1.e-8) {
                  // the grad points along MD
                  double scale = - distance_to_plane
                         /(MDnorm*sqrt(r_intersect*r_intersect
                                       + distance_to_plane*distance_to_plane));
                  *grad = MD[side] * scale;
                }
              else {
                  SCVector3 point_on_ring;
                  point_on_ring = XM_in_plane
                                * (r_intersect/XM_in_plane_norm) + M;
                  SCVector3 gradv = Xv - point_on_ring;
                  gradv.normalize();
                  *grad = gradv;
                }
            }
          distance_to_ring_in_plane =
                         r_intersect - sqrt(XM_in_plane.dot(XM_in_plane));
          return sqrt(distance_to_ring_in_plane*distance_to_ring_in_plane
                      +distance_to_plane*distance_to_plane);
        }
    }

  if (verbose) ExEnv::outn() << "returning -1.0" << endl;
  return -1.0;
}

void
Uncapped5SphereExclusionShape::boundingbox(double valuemin, double valuemax,
                                           SCVector3& p1,
                                           SCVector3& p2)
{
  SCVector3 p11;
  SCVector3 p12;
  SCVector3 p21;
  SCVector3 p22;
  SCVector3 p31;
  SCVector3 p32;

  _s1.boundingbox(valuemin,valuemax,p11,p12);
  _s2.boundingbox(valuemin,valuemax,p21,p22);
  _s3.boundingbox(valuemin,valuemax,p31,p32);

  int i;
  for (i=0; i<3; i++) {
      if (p11[i] < p21[i]) p1[i] = p11[i];
      else p1[i] = p21[i];
      if (p31[i] < p1[i]) p1[i] = p31[i];
      if (p12[i] > p22[i]) p2[i] = p12[i];
      else p2[i] = p22[i];
      if (p32[i] > p2[i]) p2[i] = p32[i];
    }
}

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

/////////////////////////////////////////////////////////////////////
// Unionshape

static ClassDesc UnionShape_cd(
  typeid(UnionShape),"UnionShape",1,"public Shape",
  0, 0, 0);

UnionShape::UnionShape()
{
}

UnionShape::~UnionShape()
{
}

void
UnionShape::add_shape(Ref<Shape> s)
{
  _shapes.insert(s);
}

// NOTE: this underestimates the distance to the surface when
//inside the surface
double
UnionShape::distance_to_surface(const SCVector3&p,SCVector3* grad) const
{
  std::set<Ref<Shape> >::const_iterator imin = _shapes.begin();
  if (imin == _shapes.end()) return 0.0;
  double min = (*imin)->distance_to_surface(p);
  for (std::set<Ref<Shape> >::const_iterator i=imin; i!=_shapes.end(); i++) {
      double d = (*i)->distance_to_surface(p);
      if (min <= 0.0) {
          if (d < 0.0 && d > min) { min = d; imin = i; }
        }
      else {
          if (min > d) { min = d; imin = i; }
        }
    }

  if (grad) {
      (*imin)->distance_to_surface(p,grad);
    }
  return min;
}

int
UnionShape::is_outside(const SCVector3&p) const
{
  for (std::set<Ref<Shape> >::const_iterator i=_shapes.begin();
       i!=_shapes.end(); i++) {
      if (!(*i)->is_outside(p)) return 0;
    }

  return 1;
}

void
UnionShape::boundingbox(double valuemin, double valuemax,
                        SCVector3& p1,
                        SCVector3& p2)
{
  if (_shapes.begin() == _shapes.end()) {
      for (int i=0; i<3; i++) p1[i] = p2[i] = 0.0;
      return;
    }
  
  SCVector3 pt1;
  SCVector3 pt2;
  
  std::set<Ref<Shape> >::iterator j = _shapes.begin();
  int i;
  (*j)->boundingbox(valuemin,valuemax,p1,p2);
  for (j++; j!=_shapes.end(); j++) {
      (*j)->boundingbox(valuemin,valuemax,pt1,pt2);
      for (i=0; i<3; i++) {
          if (pt1[i] < p1[i]) p1[i] = pt1[i];
          if (pt2[i] > p2[i]) p2[i] = pt2[i];
        }
    }
}

int
UnionShape::gradient_implemented() const
{
  for (std::set<Ref<Shape> >::const_iterator j=_shapes.begin();
       j!=_shapes.end(); j++) {
      if (!(*j)->gradient_implemented()) return 0;
    }
  return 1;
}

/////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: