molshape.cc

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

      return;
    }

  double r = sphere[0].radius() - probe_r;
  SCVector3 v1(sphere[0].x() - r, sphere[0].y() - r, sphere[0].z() - r);
  SCVector3 v2(sphere[0].x() + r, sphere[0].y() + r, sphere[0].z() + r);

  for (i=1; i<n_spheres; i++) {
      double r = sphere[i].radius() - probe_r;
      for (j=0; j<3; j++) {
          if (v1[j] > sphere[i].center()[j] - r) {
              v1[j] = sphere[i].center()[j] - r;
            }
          if (v2[j] < sphere[i].center()[j] + r) {
              v2[j] = sphere[i].center()[j] + r;
            }
        }
    }
  
  for (i=0; i<3; i++) {
      p1[i] = v1[i] - 0.01;
      p2[i] = v2[i] + 0.01;
    }
}

////////////////////////////////////////////////////////////////////////
// interval class needed by CS2Sphere

// Simple class to keep track of regions along an interval
class interval
{
    int _nsegs;                // # disjoint segments in interval
    int _max_segs;              // # segments currently allocated

    double *_min, *_max;        // arrays of ranges for segments

  private:
    // internal member function to compact interval list--this 
    // assumes that new segment is located in last element of 
    // _min and _max
    void compact(void)
    {

        if (_nsegs==1) return;

        // case 0 new segment is disjoint and below all other segments
        if (_max[_nsegs-1] < _min[0]) 
        {
            double mintmp=_min[_nsegs-1];
            double maxtmp=_max[_nsegs-1];
            for (int i=_nsegs-2; i>=0 ; i--)
            {
                _min[i+1]=_min[i];
                _max[i+1]=_max[i];
            }
            _min[0]=mintmp;
            _max[0]=maxtmp;
            return;
        }

        // case 1: new segment is disjoint and above all other segments
        if (_min[_nsegs-1] > _max[_nsegs-2]) return;

        // Fast forward to where this interval belongs
        int icount=0;
        while (_min[_nsegs-1] > _max[icount]) icount++;

        // case 2: new segment is disjoint and between two segments
        if (_max[_nsegs-1] < _min[icount]) 
        {
            double mintmp=_min[_nsegs-1];
            double maxtmp=_max[_nsegs-1];
            for (int i=_nsegs-2; i >= icount; i--)
            {
                _min[i+1]=_min[i];
                _max[i+1]=_max[i];
            }
            _min[icount]=mintmp;
            _max[icount]=maxtmp;
            return;
        }

        // new segment must overlap lower part of segment icount,
        // so redefine icount's lower boundary
        _min[icount] = (_min[_nsegs-1] < _min[icount])?
            _min[_nsegs-1]:_min[icount];

        // Now figure how far up this new segment extends
        // case 3: if it doesn't extend beyond this segment, just exit
        if (_max[_nsegs-1] < _max[icount]) { _nsegs--; return;}

        // Search forward till we find its end
        int jcount=icount;
        while (_max[_nsegs-1] > _max[jcount]) jcount++;
        
        // Case 4
        // The new segment goes to the end of all the other segments
        if (jcount == _nsegs-1)
        {
            _max[icount]=_max[_nsegs-1];
            _nsegs=icount+1;
            return;
        }

        // Case 5 
        // The new segment ends between segments
        if (_max[_nsegs-1] < _min[jcount])
        {
            _max[icount]=_max[_nsegs-1];
            // Now clobber all the segments covered by the new one
            int kcount=icount+1;
            for (int i=jcount; i<_nsegs; i++)
            {
                _min[kcount]=_min[i];
                _max[kcount]=_max[i];
                kcount++;
            }
            _nsegs=kcount-1;
            return;
        }
        
        // Case 6 
        // The new segment ends inside a segment
        if (_max[_nsegs-1] >= _min[jcount])
        {
            _max[icount]=_max[jcount];
            // Now clobber all the segments covered by the new one
            int kcount=icount+1;
            for (int i=jcount+1; i<_nsegs; i++)
            {
                _min[kcount]=_min[i];
                _max[kcount]=_max[i];
                kcount++;
            }
            _nsegs=kcount-1;
            return;
        }

        // Shouldn't get here!
        ExEnv::err0() << indent
             << "Found no matching cases in interval::compact()\n";
        print();
        exit(1);
    }

  public:
    interval(void):_nsegs(0),_max_segs(10) 
   { _min = (double*) malloc(_max_segs*sizeof(double));   // Use malloc so
      _max = (double*) malloc(_max_segs*sizeof(double));} //we can use realloc

    ~interval() { free(_min); free(_max); }
    
    // add a new segment to interval
    void add(double min, double max)
    {
        if (min > max) {double tmp=min; min=max; max=tmp;}
        if (_nsegs == _max_segs)
        {
            _max_segs *= 2;
            _min=(double *)realloc(_min, _max_segs*sizeof(double));
            _max=(double *)realloc(_max, _max_segs*sizeof(double));
        }
        
        _min[_nsegs]=min;
        _max[_nsegs]=max;
        _nsegs++;
        compact();
    }
    
    // Test to see if the interval is complete over {min, max}
    int test_interval(double min, double max)
    {
        if (_nsegs == 0) return 0;
  
        if (min > max) {double tmp=min; min=max; max=tmp;}
        
        if (min < _min[0] || max > _max[_nsegs-1]) return 0;
        for (int i=0; i < _nsegs; i++)
        {
            if (min > _min[i] && max < _max[i]) return 1;
            if (max < _min[i]) return 0;
        }
        return 0;
    }
    
    // Print out the currect state of the interval
    void print()
    {
        ExEnv::out0() << indent
             << scprintf(" _nsegs=%d; _max_segs=%d\n",_nsegs, _max_segs);
        for (int i=0; i<_nsegs; i++)
            ExEnv::out0() << indent
                 << scprintf("min[%d]=%7.4lf, max[%d]=%7.4lf\n",
                             i,_min[i],i,_max[i]); 
    }

    void clear() { _nsegs = 0; }
};

////////////////////////////////////////////////////////////////////////
// CS2Sphere

#if COUNT_CONNOLLY
int CS2Sphere::n_no_spheres_ = 0;
int CS2Sphere::n_probe_enclosed_by_a_sphere_ = 0;
int CS2Sphere::n_probe_center_not_enclosed_ = 0;
int CS2Sphere::n_surface_of_s0_not_covered_ = 0;
int CS2Sphere::n_plane_totally_covered_ = 0;
int CS2Sphere::n_internal_edge_not_covered_ = 0;
int CS2Sphere::n_totally_covered_ = 0;
#endif

void
CS2Sphere::print_counts(ostream& os)
{
  os << indent << "CS2Sphere::print_counts():\n" << incindent;
#if COUNT_CONNOLLY
  os
     << indent << "n_no_spheres = " << n_no_spheres_ << endl
     << indent << "n_probe_enclosed_by_a_sphere = "
               << n_probe_enclosed_by_a_sphere_ << endl
     << indent << "n_probe_center_not_enclosed = "
               << n_probe_center_not_enclosed_ << endl
     << indent << "n_surface_of_s0_not_covered = "
               << n_surface_of_s0_not_covered_ << endl
     << indent << "n_plane_totally_covered_ = "
               << n_plane_totally_covered_ << endl
     << indent << "n_internal_edge_not_covered = "
               << n_internal_edge_not_covered_ << endl
     << indent << "n_totally_covered = " << n_totally_covered_ << endl
     << decindent;
#else
  os << indent << "No count information is available.\n"
     << decindent;
#endif
}

// Function to determine if the centers of a bunch of spheres are separated
// by a plane from the center of another plane

// s0 is assumed to be at the origin.

// Return 1 if all of the points can be placed on the same side of a
// plane passing through s0's center.
static int
same_side(const CS2Sphere& s0, CS2Sphere *s, int n_spheres)
{
  if (n_spheres <= 3) return 1;

  SCVector3 perp;
  int sign;

  for (int i=0; i<n_spheres; i++)
    {
      for (int j=0; j<i; j++)
        {
          perp = s[i].center().perp_unit(s[j].center());
          int old_sign=0;
          for (int k=0; k < n_spheres; k++)
            {
              if (i != k && j != k)
                {
                  sign=(perp.dot(s[k].center()) < 0)? -1:1;
                  if (old_sign && old_sign != sign)
                      goto next_plane;
                  old_sign=sign;
                }
            }
          // We found a  plane with all centers on one side
          return 1;
          next_plane:
          continue;
        }
    }
  // All of the planes had points on both sides.
  return 0;
}

double
CS2Sphere::common_radius(CS2Sphere &asphere)
{
  double d=distance(asphere);
  double s=0.5*(d+_radius+asphere._radius);
  double p = s*(s-d)*(s-_radius)*(s-asphere._radius);
  //printf("common_radius: p = %5.3f\n", p);
  if (p <= 0.0) return 0.0;
  return 2.*sqrt(p)/d;
}

#define PRINT_SPECIAL_CASES 0
#if PRINT_SPECIAL_CASES
static void
print_spheres(const CS2Sphere& s0, CS2Sphere* s, int n_spheres)
{
  static int output_number;
  char filename[80];
  sprintf(filename,"spherelist_%d.oogl",output_number);
  FILE* fp = fopen(filename,"w");
  fprintf(fp,"LIST\n");
  fprintf(fp,"{\n");
  fprintf(fp,"  appearance {\n");
  fprintf(fp,"      material {\n");
  fprintf(fp,"         ambient 0.5 0.1 0.1\n");
  fprintf(fp,"         diffuse 1.0 0.2 0.2\n");
  fprintf(fp,"       }\n");
  fprintf(fp,"    }\n");
  fprintf(fp," = SPHERE\n");
  fprintf(fp," %15.8f %15.8f %15.8f %15.8f\n",
          s0.radius(), s0.x(), s0.y(), s0.z());
  fprintf(fp,"}\n");
  for (int i=0; i<n_spheres; i++) {
      fprintf(fp,"{ = SPHERE\n");
      fprintf(fp," %15.8f %15.8f %15.8f %15.8f\n",
              s[i].radius(), s[i].x(), s[i].y(), s[i].z());
      fprintf(fp,"}\n");
    }
  fclose(fp);
  output_number++;
}
#endif

// Function to determine if there is any portion of s0 that 
// is not inside one or more of the spheres in s[]
int
CS2Sphere::intersect(CS2Sphere *s, int n_spheres) const
{
    if (n_spheres == 0) {
        n_no_spheres_++;
        return 0;
      }
    CS2Sphere s0;
    s0 = *this;
    // Declare an interval object to manage overlap information
    // it is static so it will only call malloc twice
    static interval intvl;
    // First make sure that at least one sphere in s[] contains
    // the center of s0 and that s0 is not contained inside
    // one of the spheres
    int center_is_contained = 0;
    int i;
    for (i=0; i<n_spheres; i++)
    {
        double d=s0.distance(s[i]);
        if (d+s0.radius() < s[i].radius()) {
            n_probe_enclosed_by_a_sphere_++;
            return 1;
          }
        if (d < s[i].radius()) center_is_contained = 1;
    }
    if (!center_is_contained) {
        n_probe_center_not_enclosed_++;
        return 0;
      }
    
    // Let's first put s0 at the origin
    for (i=0; i<n_spheres; i++)
        s[i].recenter(s0.center());
    s0.recenter(s0.center());
            
    // Now check to make sure that the surface of s0 is completely
    // included in spheres in s[], by making sure that all the
    // circles describing the intersections of every sphere with
    // s0 are included in at least one other sphere.
    double epsilon=1.e-8;
    for (i=0; i<n_spheres; i++)
    {
        // calculate radius of the intersection of s0 and s[i]
        double cr = s0.common_radius(s[i]);
        if (cr == 0.0) {
            continue;
        }
        
        // We're chosing that the intersection of s[i] and s0 
        // occurs parallel to the x-y plane, so we'll need to rotate the
        // center of s[j] appropriately.  
        // Create a rotation matrix that take the vector from 
        // the centers of s0 to s[i] and puts it on the z axis
        static const SCVector3 Zaxis(0.0, 0.0, 1.0);
        SCMatrix3 rot = rotation_mat(s0.center_vec(s[i]),Zaxis);
        
        // Now calculate the Z position of the intersection of 
        // s0 and s[i]
        double d=s0.distance(s[i]);
        double z_plane;
        if (s[i].radius()*s[i].radius() < d*d+s0.radius()*s0.radius())
            z_plane=sqrt(s0.radius()*s0.radius()-cr*cr);
        else
            z_plane=-sqrt(s0.radius()*s0.radius()-cr*cr);
        
        // Initialize the interval object
        intvl.clear();
        
        // Loop over the other spheres
        for (int j=0; j<n_spheres; j++)
            if (i != j)
            {
                // Rotate the center of s[j] to appropriate refence frame
                SCVector3 rcent = rot*s0.center_vec(s[j]);
                
                double x0=rcent.x();
                double y0=rcent.y();
                double z0=rcent.z();
                
                // Does this sphere even reach the plane where
                // the intersection of s0 and s[i] occurs?
                // If not, let's go to the next sphere