shape.cc

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

  o << decindent;
  o << indent << "s2 = ";
  o << incindent << skipnextindent;
  _s2.print(o);
  o << decindent;
  o << decindent;
}

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

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

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

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

/////////////////////////////////////////////////////////////////////
// is in triangle

static int
is_in_unbounded_triangle(const SCVector3&XP,const SCVector3&AP,const SCVector3&BP)
{
  SCVector3 axis = BP.cross(AP);

  SCVector3 BP_perp = BP; BP_perp.rotate(M_PI_2,axis);
  double u = BP_perp.dot(XP)/BP_perp.dot(AP);
  if (u<0.0) return 0;

  SCVector3 AP_perp = AP; AP_perp.rotate(M_PI_2,axis);
  double w = AP_perp.dot(XP)/AP_perp.dot(BP);
  if (w<0.0) return 0;

  return 1;
}

/////////////////////////////////////////////////////////////////////
// ReentrantUncappedTorusHoleShape

static ClassDesc ReentrantUncappedTorusHoleShape_cd(
  typeid(ReentrantUncappedTorusHoleShape),"ReentrantUncappedTorusHoleShape",1,"public UncappedTorusHoleShape",
  0, 0, 0);

ReentrantUncappedTorusHoleShape::ReentrantUncappedTorusHoleShape(double r,
                                                 const SphereShape& s1,
                                                 const SphereShape& s2):
  UncappedTorusHoleShape(r,s1,s2)
{
  rAP = r + s1.radius();
  rBP = r + s2.radius();
  BA = B() - A();

  // Find the points at the ends of the two cone-like objects, I[0] and I[1].
  // they are given by:
  //   I = z BA, where BA = B-A and I is actually IA = I - A
  //   r^2 = PI.PI, where PI = PA-I and P is the center of a probe sphere
  // this gives
  //  z^2 BA.BA - 2z PA.BA + PA.PA - r^2 = 0

  SCVector3 arbitrary; 
  arbitrary[0] = arbitrary[1] = arbitrary[2] = 0.0;
  SCVector3 P;
  in_plane_sphere(arbitrary,P);
  SCVector3 PA = P - A();

  double a = BA.dot(BA);
  double minus_b = 2.0 * PA.dot(BA);
  double c = PA.dot(PA) - r*r;
  double b2m4ac = minus_b*minus_b - 4*a*c;
  double sb2m4ac;
  if (b2m4ac >= 0.0) {
      sb2m4ac = sqrt(b2m4ac);
    }
  else if (b2m4ac > -1.0e-10) {
      sb2m4ac = 0.0;
    }
  else {
      ExEnv::errn() << "ReentrantUncappedTorusHoleShape:: imaginary point" << endl;
      abort();
    }
  double zA = (minus_b - sb2m4ac)/(2.0*a);
  double zB = (minus_b + sb2m4ac)/(2.0*a);
  I[0] = BA * zA + A();
  I[1] = BA * zB + A();
}
ReentrantUncappedTorusHoleShape::~ReentrantUncappedTorusHoleShape()
{
}
int
ReentrantUncappedTorusHoleShape::
  is_outside(const SCVector3&X) const
{
  SCVector3 Xv(X);

  SCVector3 P;
  in_plane_sphere(X,P);
  SCVector3 XP = Xv - P;
  double rXP = XP.norm();
  if (rXP > rAP || rXP > rBP) return 1;

  SCVector3 AP = A() - P;
  SCVector3 BP = B() - P;

  if (!is_in_unbounded_triangle(XP,AP,BP)) return 1;

  if (rXP < radius()) {
      return 1;
    }

  return 0;
}
double
ReentrantUncappedTorusHoleShape::
  distance_to_surface(const SCVector3&X,SCVector3*grad) const
{
  SCVector3 Xv(X);

  SCVector3 P;
  in_plane_sphere(X,P);
  SCVector3 XP = Xv - P;
  double rXP = XP.norm();
  if (rXP > rAP || rXP > rBP) return shape_infinity;

  SCVector3 AP = A() - P;
  SCVector3 BP = B() - P;

  if (!is_in_unbounded_triangle(XP,AP,BP)) return shape_infinity;

  SCVector3 I1P = I[0] - P;
  SCVector3 I2P = I[1] - P;

  if (!is_in_unbounded_triangle(XP,I1P,I2P)) {
      if (rXP < radius()) {
          if (grad) {
              SCVector3 unit(XP);
              unit.normalize();
              *grad = unit;
            }
          return radius() - rXP;
        }
      else return -1.0;
    }

  return closest_distance(Xv,(SCVector3*)I,2,grad);
}

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

/////////////////////////////////////////////////////////////////////
// NonreentrantUncappedTorusHoleShape

static ClassDesc NonreentrantUncappedTorusHoleShape_cd(
  typeid(NonreentrantUncappedTorusHoleShape),"NonreentrantUncappedTorusHoleShape",1,"public UncappedTorusHoleShape",
  0, 0, 0);

NonreentrantUncappedTorusHoleShape::
  NonreentrantUncappedTorusHoleShape(double r,
                                     const SphereShape& s1,
                                     const SphereShape& s2):
  UncappedTorusHoleShape(r,s1,s2)
{
  rAP = r + s1.radius();
  rBP = r + s2.radius();
  BA = B() - A();
}
NonreentrantUncappedTorusHoleShape::~NonreentrantUncappedTorusHoleShape()
{
}
double NonreentrantUncappedTorusHoleShape::
  distance_to_surface(const SCVector3&X,SCVector3* grad) const
{
  SCVector3 Xv(X);

  SCVector3 P;
  in_plane_sphere(X,P);
  SCVector3 PX = P - Xv;
  double rPX = PX.norm();
  if (rPX > rAP || rPX > rBP) return shape_infinity;

  SCVector3 PA = P - A();
  SCVector3 XA = Xv - A();

  SCVector3 axis = BA.cross(PA);

  SCVector3 BA_perp = BA; BA_perp.rotate(M_PI_2,axis);
  double u = BA_perp.dot(XA)/BA_perp.dot(PA);
  if (u<0.0 || u>1.0) return shape_infinity;

  SCVector3 PA_perp = PA; PA_perp.rotate(M_PI_2,axis);
  double w = PA_perp.dot(XA)/PA_perp.dot(BA);
  if (w<0.0 || w>1.0) return shape_infinity;

  double uw = u+w;
  if (uw<0.0 || uw>1.0) return shape_infinity;

  if (rPX < radius()) {
      if (grad) {
          SCVector3 unit(PX);
          unit.normalize();
          *grad = unit;
        }
      return radius() - rPX;
    }

  return -1;
}

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

/////////////////////////////////////////////////////////////////////
// Uncapped5SphereExclusionShape

static ClassDesc Uncapped5SphereExclusionShape_cd(
  typeid(Uncapped5SphereExclusionShape),"Uncapped5SphereExclusionShape",1,"public Shape",
  0, 0, 0);

Uncapped5SphereExclusionShape*
Uncapped5SphereExclusionShape::
  newUncapped5SphereExclusionShape(double r,
                                   const SphereShape& s1,
                                   const SphereShape& s2,
                                   const SphereShape& s3)
{
  Uncapped5SphereExclusionShape* s =
    new Uncapped5SphereExclusionShape(r,s1,s2,s3);
  if (s->solution_exists()) {
      return s;
    }
  delete s;
  return 0;
}
static int verbose = 0;
Uncapped5SphereExclusionShape::
  Uncapped5SphereExclusionShape(double radius,
                                const SphereShape&s1,
                                const SphereShape&s2,
                                const SphereShape&s3):
  _s1(s1),
  _s2(s2),
  _s3(s3),
  _r(radius)
{
  double rAr = rA() + r();
  double rAr2 = rAr*rAr;
  double rBr = rB() + r();
  double rBr2 = rBr*rBr;
  double rCr = rC() + r();
  double rCr2 = rCr*rCr;
  double A2 = A().dot(A());
  double B2 = B().dot(B());
  double C2 = C().dot(C());
  SCVector3 BA = B()-A();
  double DdotBA = 0.5*(rAr2 - rBr2 + B2 - A2);
  double DAdotBA = DdotBA - A().dot(BA);
  double BA2 = BA.dot(BA);
  SCVector3 CA = C() - A();
  double CAdotBA = CA.dot(BA);
  SCVector3 CA_perpBA = CA - (CAdotBA/BA2)*BA;
  double CA_perpBA2 = CA_perpBA.dot(CA_perpBA);
  if (CA_perpBA2 < 1.0e-23) {
      _solution_exists = 0;
      return;
    }
  double DdotCA_perpBA = 0.5*(rAr2 - rCr2 + C2 - A2)
    - CAdotBA*DdotBA/BA2;
  double DAdotCA_perpBA = DdotCA_perpBA - A().dot(CA_perpBA);
  double rAt2 = DAdotBA*DAdotBA/BA2 + DAdotCA_perpBA*DAdotCA_perpBA/CA_perpBA2;
  double h2 = rAr2 - rAt2;
  if (h2 <= 0.0) {
      _solution_exists = 0;
      return;
    }
  _solution_exists = 1;

  double h = sqrt(h2);
  if (h<r()) {
      _reentrant = 1;
      //ExEnv::outn() << "WARNING: throwing out reentrant shape" << endl;
      //_solution_exists = 0;
      //return;
    }
  else {
      _reentrant = 0;
      //ExEnv::outn() << "WARNING: throwing out nonreentrant shape" << endl;
      //_solution_exists = 0;
      //return;
    }

  // The projection of D into the ABC plane
  SCVector3 MA = (DAdotBA/BA2)*BA + (DAdotCA_perpBA/CA_perpBA2)*CA_perpBA;
  M = MA + A();
  SCVector3 BAxCA = BA.cross(CA);
  D[0] = M + h * BAxCA * ( 1.0/BAxCA.norm() );
  D[1] = M - h * BAxCA * ( 1.0/BAxCA.norm() );

  // The projection of D into the ABC plane, M, does not lie in the
  // ABC triangle, then this shape must be treated carefully and it
  // will be designated as folded.
  SCVector3 MC = M - C();
  if (!(is_in_unbounded_triangle(MA, BA, CA)
        &&is_in_unbounded_triangle(MC, B() - C(), A() - C()))) {
      _folded = 1;
      SCVector3 MB = M - B();
      double MA2 = MA.dot(MA);
      double MB2 = MB.dot(MB);
      double MC2 = MC.dot(MC);
      if (MA2 < MB2) {
          F1 = A();
          if (MB2 < MC2) F2 = B();
          else F2 = C();
        }
      else {
          F1 = B();
          if (MA2 < MC2) F2 = A();
          else F2 = C();
        }
    }
  else _folded = 0;
  
  //ExEnv::outn() << scprintf("r = %14.8f, h = %14.8f\n",r(),h);
  //M.print();
  //D[0].print();
  //D[1].print();
  //A().print();
  //B().print();
  //C().print();

  int i;
  for (i=0; i<2; i++) {
      SCVector3 AD = A() - D[i];
      SCVector3 BD = B() - D[i];
      SCVector3 CD = C() - D[i];
      BDxCD[i] = BD.cross(CD);
      BDxCDdotAD[i] = BDxCD[i].dot(AD);
      CDxAD[i] = CD.cross(AD);
      CDxADdotBD[i] = CDxAD[i].dot(BD);
      ADxBD[i] = AD.cross(BD);
      ADxBDdotCD[i] = ADxBD[i].dot(CD);
    }

  for (i=0; i<2; i++) MD[i] = M - D[i];

  // reentrant surfaces need a whole bunch more to be able to compute
  // the distance to the surface
  if (_reentrant) {
      int i;
      double rMD = MD[0].norm(); // this is the same as rMD[1]
      theta_intersect = M_PI_2 - asin(rMD/r());
      r_intersect = r() * sin(theta_intersect);

      {
        // Does the circle of intersection intersect with BA?
        SCVector3 MA = M - A();
        SCVector3 uBA = B() - A(); uBA.normalize();
        SCVector3 XA = uBA * MA.dot(uBA);
        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_AB = 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++) {
                    IABD[i][j] = XM + d_intersect_from_x[j]*uBA + MD[i];