molshape.cc

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

                double z_dist=s[j].radius()*s[j].radius()-
                    (z0-z_plane)*(z0-z_plane);
                if (z_dist < 0.0)
                    continue;
                
                // Calculate radius of circular projection of s[j]
                // onto s0-s[i] intersection plane
                double r_2=z_dist;
                
                // Precalculate a bunch of factors 
                double cr_2=cr*cr;
                double x0_2=x0*x0; double y0_2=y0*y0;
                double dist=sqrt(x0_2+y0_2);
                
                // If the projection of s[j] on x-y doesn't reach the
                // intersection of s[i] and s0, continue.
                if (r_2 < (dist-cr)*(dist-cr))
                    continue;
                
                // If the projection of s[j] on x-y engulfs the intersection
                // of s[i] and s0, cover interval and continue
                if (r_2 > (dist+cr)*(dist+cr))
                {
                    intvl.add(0, 2.*M_PI);
                    continue;
                }
                
                // Calculation the radical in the quadratic equation
                // determining the overlap of the two circles
                double radical=x0_2*(-cr_2*cr_2 + 2*cr_2*r_2 - 
                                     r_2*r_2 + 2*cr_2*x0_2 + 
                                     2*r_2*x0_2 - x0_2*x0_2 + 
                                     2*cr_2*y0_2 + 2*r_2*y0_2 - 
                                     2*x0_2*y0_2 - y0_2*y0_2);
                
                // Check to see if there's any intersection at all
                // I.e. if one circle is inside the other  (Note that
                // we've already checked to see if s[j] engulfs
                // the intersection of s0 and s[i])
                if (radical <= 0.0) continue;
                
                // Okay, go ahead and calculate the intersection points
                double x_numer = cr_2*x0_2 - r_2*x0_2 + x0_2*x0_2 + x0_2*y0_2;
                double x_denom = 2*x0*x0_2 + 2*x0*y0_2;
                double y_numer = cr_2*y0 - r_2*y0 + x0_2*y0 + y0*y0_2;
                double y_denom = 2*(x0_2 + y0_2);
                
                double sqrt_radical = sqrt(radical);
                
                double x_0=(x_numer - y0*sqrt_radical)/x_denom;
                double y_0=(y_numer + sqrt_radical)/y_denom;
                double x_1=(x_numer + y0*sqrt_radical)/x_denom;
                double y_1=(y_numer - sqrt_radical)/y_denom;
                
                // Now calculate the angular range of these ordered
                // points and place them on the first Riemann sheet.
                // and sort their order
                double theta1=atan2(y_0, x_0);
                double theta2=atan2(y_1, x_1);
                if (theta1 < 0.0) theta1+=2.*M_PI;
                if (theta2 < 0.0) theta2+=2.*M_PI;
                if (theta1 > theta2)
                {
                    double tmptheta=theta1;
                    theta1=theta2;
                    theta2=tmptheta;
                }
                
                // Determine which of the two possible chords 
                // is inside s[j]
                double dor=(x0-cr)*(x0-cr)+y0*y0;
                if (dor < r_2)
                {
                    intvl.add(0, theta1);
                    intvl.add(theta2, 2.*M_PI);
                }
                else            
                {
                    intvl.add(theta1, theta2);
                }
                
                // Now test to see if the range is covered
                if (intvl.test_interval(epsilon, 2.*M_PI-epsilon))
                {
                    // No need to keep testing, move on to next i
                    break;
                }
                
            }
        // If the intersection wasn't totally covered, the sphere
        // intersection is incomplete
        if (!intvl.test_interval(epsilon, 2.*M_PI-epsilon)) {
            n_surface_of_s0_not_covered_++;
            // goto next_test;
            return 0;
        }
    }

    // for the special case of all sphere's centers on one side of
    // a plane passing through s0's center we are done; the probe
    // must be completely intersected.
    if (same_side(s0,s,n_spheres)) {
        n_plane_totally_covered_++;
        return 1;
      }
    
    // As a final test of the surface coverage, make sure that all
    // of the intersection surfaces between s0 and s[] are included 
    // inside more than one sphere.
    int angle_segs;
    double max_angle[2], min_angle[2];
    for (i=0; i<n_spheres; i++)
    {
        // For my own sanity, let's put s[i] at the origin 
        int k;
        for (k=0; k<n_spheres; k++)
            if (k != i)
                s[k].recenter(s[i].center());
        s0.recenter(s[i].center());
        s[i].recenter(s[i].center());
        
        for (int j=0; j<i; j++)
        {

            // calculate radius of the intersection of s[i] and s[j]
            double cr = s[i].common_radius(s[j]);
            if (cr == 0.0) {
                continue;                   // s[i] and s[j] don't intersect
            }

            // We're chosing that the intersection of s[i] and s[j]
            // occurs parallel to the x-y plane, so we'll need to rotate the
            // center of all s[]'s and s0 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(s[i].center_vec(s[j]),Zaxis);
            
            // Now calculate the Z position of the intersection of 
            // s[i] and s[j]
            double d=s[i].distance(s[j]);
            double z_plane;
            if (s[j].radius()*s[j].radius() < s[i].radius()*s[i].radius()+d*d)
                z_plane=sqrt(s[i].radius()*s[i].radius()-cr*cr);
            else
                z_plane=-sqrt(s[i].radius()*s[i].radius()-cr*cr);
            
            // Determine which part of the this intersection
            // occurs within s0
            // Rotate the center of s0 to appropriate refence frame
            SCVector3 rcent = rot*s[i].center_vec(s0);
            
            double x0=rcent.x();
            double y0=rcent.y();
            double z0=rcent.z();
            
            // Does this s0 even reach the plane where
            // the intersection of s[i] and s[j] occurs?
            // If not, let's go to the next sphere j
            double z_dist=s0.radius()*s0.radius()-
                (z0-z_plane)*(z0-z_plane);
            if (z_dist < 0.0)
                continue;
            
            // Calculate radius of circular projection of s0
            // onto s[i]-s[j] intersection plane
            double r_2=z_dist;
            
            // Precalculate a bunch of factors 
            double cr_2=cr*cr;
            double x0_2=x0*x0; double y0_2=y0*y0;
            double dist=sqrt(x0_2+y0_2);
            
            // If the projection of s[j] on x-y doesn't reach the
            // intersection of s[i] and s0, continue.
            if (r_2 < (dist-cr)*(dist-cr))
                continue;
            
            // If the projection of s0 on x-y engulfs the intersection
            // of s[i] and s[j], the intersection interval is 0 to 2pi
            if (r_2 > (dist+cr)*(dist+cr))
            {
                angle_segs=1;
                min_angle[0]=0.0;
                max_angle[0]=2.*M_PI;
            }
            
            // Calculation the radical in the quadratic equation
            // determining the overlap of the two circles
            double radical=x0_2*(-cr_2*cr_2 + 2*cr_2*r_2 - 
                                 r_2*r_2 + 2*cr_2*x0_2 + 
                                 2*r_2*x0_2 - x0_2*x0_2 + 
                                 2*cr_2*y0_2 + 2*r_2*y0_2 - 
                                 2*x0_2*y0_2 - y0_2*y0_2);
            
            // Check to see if there's any intersection at all
            // I.e. if one circle is inside the other  (Note that
            // we've already checked to see if s0 engulfs
            // the intersection of s[i] and s[j]), so this
            // must mean that the intersection of s[i] and s[j]
            // occurs outside s0
            if (radical <= 0.0) continue;

            // Okay, go ahead and calculate the intersection points
            double x_numer = cr_2*x0_2 - r_2*x0_2 + x0_2*x0_2 + x0_2*y0_2;
            double x_denom = 2*x0*x0_2 + 2*x0*y0_2;
            double y_numer = cr_2*y0 - r_2*y0 + x0_2*y0 + y0*y0_2;
            double y_denom = 2*(x0_2 + y0_2);
            
            double sqrt_radical = sqrt(radical);
            
            double x_0=(x_numer - y0*sqrt_radical)/x_denom;
            double y_0=(y_numer + sqrt_radical)/y_denom;
            double x_1=(x_numer + y0*sqrt_radical)/x_denom;
            double y_1=(y_numer - sqrt_radical)/y_denom;
            
            // Now calculate the angular range of these ordered
            // points and place them on the first Riemann sheet.
            // and sort their order
            double theta1=atan2(y_0, x_0);
            double theta2=atan2(y_1, x_1);
            if (theta1 < 0.0) theta1+=2.*M_PI;
            if (theta2 < 0.0) theta2+=2.*M_PI;
            if (theta1 > theta2)
            {
                double tmptheta=theta1;
                theta1=theta2;
                theta2=tmptheta;
            }
            //printf("theta1=%lf, theta2=%lf\n",theta1,theta2);

            // Determine which of the two possible chords 
            // is inside s0

            // But first see if s0 is inside this intersection:
            double origin_dist=((x0-cr)*(x0-cr)+(y0*y0));
            if (origin_dist < r_2) // it's the angle containing
                                       // the origin
            {
                angle_segs=2;
                min_angle[0]=0.0;
                max_angle[0]=theta1;
                min_angle[1]=theta2;
                max_angle[1]=2.*M_PI;
            }
            else            // it's the angle not including the origin
            {
                angle_segs=1;
                min_angle[0]=theta1;
                max_angle[0]=theta2;
            }
            
            // Initialize the interval object
            intvl.clear();
            
            // Loop over the other spheres
            for (k=0; k<n_spheres; k++)
            {
                if (k != i && k != j)
                {
                    // Rotate the center of s[k] to appropriate reference frame
                    rcent = rot*s[i].center_vec(s[k]);
                    
                    double x0=rcent.x();
                    double y0=rcent.y();
                    double z0=rcent.z();
                    
                    // Does this sphere even reach the plane where
                    // the intersection of s[i] and s[j] occurs?
                    // If not, let's go to the next sphere
                    double z_dist=s[k].radius()*s[k].radius()-
                        (z0-z_plane)*(z0-z_plane);
                    if (z_dist < 0.0)
                        continue;
                    
                    // Calculate radius of circular projection of s[k]
                    // onto s[i]-s[j] intersection plane
                    double r_2=z_dist;
                    
                    // Precalculate a bunch of factors 
                    double cr_2=cr*cr;
                    double x0_2=x0*x0; double y0_2=y0*y0;
                    double dist=sqrt(x0_2+y0_2);
                    
                    // If the projection of s[k] on x-y doesn't reach the
                    // intersection of s[i] and s[j], continue.
                    if (r_2 < (dist-cr)*(dist-cr))
                        continue;
                    
                    // If the projection of s[k] on x-y engulfs the intersection
                    // of s[i] and s0, cover interval and continue
                    if (r_2 > (dist+cr)*(dist+cr))
                    {
                        intvl.add(0, 2.*M_PI);
                        continue;
                    }
                    
                    // Calculation the radical in the quadratic equation
                    // determining the overlap of the two circles
                    radical=x0_2*(-cr_2*cr_2 + 2*cr_2*r_2 - 
                                  r_2*r_2 + 2*cr_2*x0_2 + 
                                  2*r_2*x0_2 - x0_2*x0_2 + 
                                  2*cr_2*y0_2 + 2*r_2*y0_2 - 
                                  2*x0_2*y0_2 - y0_2*y0_2);
                    
                    // Check to see if there's any intersection at all
                    // I.e. if one circle is inside the other  (Note that
                    // we've already checked to see if s[k] engulfs
                    // the intersection of s[i] and s[j])
                    if (radical <= 0.0) continue;
                    
                    // Okay, go ahead and calculate the intersection points
                    x_numer = cr_2*x0_2 - r_2*x0_2 + x0_2*x0_2 + x0_2*y0_2;
                    x_denom = 2*x0*x0_2 + 2*x0*y0_2;
                    y_numer = cr_2*y0 - r_2*y0 + x0_2*y0 + y0*y0_2;
                    y_denom = 2*(x0_2 + y0_2);
                    
                    sqrt_radical = sqrt(radical);
                    
                    double x_0=(x_numer - y0*sqrt_radical)/x_denom;
                    double y_0=(y_numer + sqrt_radical)/y_denom;
                    double x_1=(x_numer + y0*sqrt_radical)/x_denom;
                    double y_1=(y_numer - sqrt_radical)/y_denom;

                    // Now calculate the angular range of these ordered
                    // points and place them on the first Riemann sheet.
                    // and sort their order
                    theta1=atan2(y_0, x_0);
                    theta2=atan2(y_1, x_1);
                    if (theta1 < 0.0) theta1+=2.*M_PI;
                    if (theta2 < 0.0) theta2+=2.*M_PI;
                    if (theta1 > theta2)
                    {
                        double tmptheta=theta1;
                        theta1=theta2;
                        theta2=tmptheta;
                    }
                    //printf("In k loop, k=%d, theta1=%lf, theta2=%lf\n",
                    //       k,theta1, theta2);
                    // Determine which of the two possible chords 
                    // is inside s[k]
                    double origin_dist=((x0-cr)*(x0-cr)+(y0*y0));
                    if (origin_dist < r_2) // it's got the origin
                    {
                        intvl.add(0, theta1);
                        intvl.add(theta2, 2.*M_PI);
                    }
                    else        // it doesn't have the origin
                    {
                        intvl.add(theta1, theta2);
                    }

                    // Now test to see if the range is covered
                    if (intvl.test_interval(min_angle[0]+epsilon,
                                            max_angle[0]-epsilon) &&
                        (angle_segs!=2 ||
                         intvl.test_interval(min_angle[1]+epsilon,
                                             max_angle[1]-epsilon)))
                    {
                        goto next_j;
                    }
                }
            }
            if (!intvl.test_interval(min_angle[0]+epsilon,
                                     max_angle[0]-epsilon))
            {
                // No need to keep testing, return 0
                n_internal_edge_not_covered_++;
                return 0;
                //printf(" Non-internal coverage(1)\n");
                //goto next_test;
            }
            if (angle_segs==2)
            {
                if  (!intvl.test_interval(min_angle[1]+epsilon,
                                          max_angle[1]-epsilon))
                {
                    n_internal_edge_not_covered_++;
                    return 0;
                    //printf(" Non-internal coverage(2)\n");
                    //goto next_test;
                }
                else
                {
                    goto next_j;
                }
            }
          next_j:
            continue;
        }
    }
    
    // Since we made it past all of the sphere intersections, the
    // surface is totally covered
    n_totally_covered_++;
    return 1;
}

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

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