molshape.cc

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

//
// molshape.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 <vector>

#include <chemistry/molecule/molshape.h>
#include <chemistry/molecule/molecule.h>
#include <math/scmat/matrix3.h>

using namespace std;
using namespace sc;

////////////////////////////////////////////////////////////////////////
// VDWShape

static ClassDesc VDWShape_cd(
  typeid(VDWShape),"VDWShape",1,"public UnionShape",
  0, create<VDWShape>, 0);

VDWShape::VDWShape(const Ref<Molecule>&mol)
{
  initialize(mol);
}

VDWShape::VDWShape(const Ref<KeyVal>&keyval)
{
  Ref<Molecule> mol; mol << keyval->describedclassvalue("molecule");
  atominfo_ << keyval->describedclassvalue("atominfo");
  initialize(mol);
}

void
VDWShape::initialize(const Ref<Molecule>&mol)
{
  Ref<AtomInfo> a;
  if (atominfo_.null()) a = mol->atominfo();
  else a = atominfo_;

  _shapes.clear();
  for (int i=0; i<mol->natom(); i++) {
      SCVector3 r;
      for (int j=0; j<3; j++) r[j] = mol->r(i,j);
      add_shape(
          new SphereShape(r,a->vdw_radius(mol->Z(i)))
          );
    }
}

VDWShape::~VDWShape()
{
}

////////////////////////////////////////////////////////////////////////
// static functions for DiscreteConnollyShape and ConnollyShape

static double
find_atom_size(const Ref<AtomInfo>& a, int Z)
{
  return a->vdw_radius(Z);
}

////////////////////////////////////////////////////////////////////////
// DiscreteConnollyShape

static ClassDesc DiscreteConnollyShape_cd(
  typeid(DiscreteConnollyShape),"DiscreteConnollyShape",1,"public UnionShape",
  0, create<DiscreteConnollyShape>, 0);

DiscreteConnollyShape::DiscreteConnollyShape(const Ref<KeyVal>&keyval)
{
  Ref<Molecule> mol; mol << keyval->describedclassvalue("molecule");
  double probe_radius = keyval->doublevalue("probe_radius");
  if (keyval->error() != KeyVal::OK) {
      probe_radius = 2.6456173;
    }
  radius_scale_factor_ = keyval->doublevalue("radius_scale_factor");
  if (keyval->error() != KeyVal::OK) {
      radius_scale_factor_ = 1.2;
    }
  atominfo_ << keyval->describedclassvalue("atominfo");
  initialize(mol,probe_radius);
}

void
DiscreteConnollyShape::initialize(const Ref<Molecule>&mol,double probe_radius)
{
  _shapes.clear();
  std::vector<Ref<SphereShape> > spheres(0);

  Ref<AtomInfo> a;
  if (atominfo_.null()) a = mol->atominfo();
  else a = atominfo_;

  int i;
  for (i=0; i<mol->natom(); i++) {
      SCVector3 r(mol->r(i));
      Ref<SphereShape>
        sphere(
            new SphereShape(r,radius_scale_factor_*find_atom_size(a,
                                             mol->Z(i)))
            );
      add_shape(sphere.pointer());
      spheres.push_back(sphere);
    }

  ////////////////////// Leave out the other shapes
  //return;

  for (i=0; i<spheres.size(); i++) {
      for (int j=0; j<i; j++) {
          Ref<Shape> th =
            UncappedTorusHoleShape::newUncappedTorusHoleShape(probe_radius,
                                              *(spheres[i].pointer()),
                                              *(spheres[j].pointer()));
          if (th.null()) continue;
          add_shape(th);

          ////////////////////// Leave out the three sphere shapes
          //continue;
          
          // now check for excluding volume for groups of three spheres
          for (int k=0; k<j; k++) {
              Ref<Shape> e =
                Uncapped5SphereExclusionShape::
              newUncapped5SphereExclusionShape(probe_radius,
                                               *(spheres[i].pointer()),
                                               *(spheres[j].pointer()),
                                               *(spheres[k].pointer()));
              if (e.nonnull()) add_shape(e);
            }
        }
    }
}

DiscreteConnollyShape::~DiscreteConnollyShape()
{
}

////////////////////////////////////////////////////////////////////////
// ConnollyShape

static ClassDesc ConnollyShape_cd(
  typeid(ConnollyShape),"ConnollyShape",1,"public Shape",
  0, create<ConnollyShape>, 0);

ConnollyShape::ConnollyShape(const Ref<KeyVal>&keyval)
{
  box_ = 0;
  sphere = 0;
  Ref<Molecule> mol; mol << keyval->describedclassvalue("molecule");
  probe_r = keyval->doublevalue("probe_radius");
  if (keyval->error() != KeyVal::OK) {
      probe_r = 2.6456173;
    }
  atominfo_ << keyval->describedclassvalue("atominfo");
  radius_scale_factor_ = keyval->doublevalue("radius_scale_factor");
  if (keyval->error() != KeyVal::OK) {
      radius_scale_factor_ = 1.2;
    }
  initialize(mol,probe_r);
}

#if COUNT_CONNOLLY
int ConnollyShape::n_total_ = 0;
int ConnollyShape::n_inside_vdw_ = 0;
int ConnollyShape::n_with_nsphere_[CONNOLLYSHAPE_N_WITH_NSPHERE_DIM];
#endif

void
ConnollyShape::print_counts(ostream& os)
{
  os << indent << "ConnollyShape::print_counts():\n" << incindent;
#if COUNT_CONNOLLY
  os
     << indent << "n_total = " << n_total_ << endl
     << indent << "n_inside_vdw = " << n_inside_vdw_ << endl;
  for (int i=0; i<CONNOLLYSHAPE_N_WITH_NSPHERE_DIM-1; i++) {
      os << indent
         << scprintf("n with nsphere = %2d: %d\n", i, n_with_nsphere_[i]);
    }
  os << indent
     << scprintf("n with nsphere >= %d: %d\n",
                 CONNOLLYSHAPE_N_WITH_NSPHERE_DIM-1,
                 n_with_nsphere_[CONNOLLYSHAPE_N_WITH_NSPHERE_DIM-1])
     << decindent;
#else
  os << indent << "No count information is available.\n" << decindent;
#endif
}

void
ConnollyShape::initialize(const Ref<Molecule>&mol,double probe_radius)
{
  clear();

  n_spheres = mol->natom();
  sphere = new CS2Sphere[n_spheres];

  Ref<AtomInfo> a;
  if (atominfo_.null()) a = mol->atominfo();
  else a = atominfo_;

  int i;
  for (i=0; i<n_spheres; i++) {
      SCVector3 r(mol->r(i));
      sphere[i].initialize(r,radius_scale_factor_*find_atom_size(a,
                                            mol->Z(i))
                              + probe_r);
    }

  // initialize a grid of lists of local spheres
  if (n_spheres) {
      // find the bounding box
      SCVector3 lower(sphere[0].center()), upper(sphere[0].center());
      for (i=0; i<n_spheres; i++) {
          SCVector3 l(sphere[i].center()), u(sphere[i].center());
          for (int j=0; j<3; j++) {
              l[j] -= probe_r + sphere[i].radius();
              u[j] += probe_r + sphere[i].radius();
              if (l[j]<lower[j]) lower[j] = l[j];
              if (u[j]>upper[j]) upper[j] = u[j];
            }
        }
      // compute the parameters for converting x, y, z into a box number
      lower_ = lower;
      l_ = 10.0;
      xmax_ = (int)((upper[0]-lower[0])/l_);
      ymax_ = (int)((upper[1]-lower[1])/l_);
      zmax_ = (int)((upper[2]-lower[2])/l_);
      // allocate the boxes
      box_ = new std::vector<int>**[xmax_+1];
      for (i=0; i<=xmax_; i++) {
          box_[i] = new std::vector<int>*[ymax_+1];
          for (int j=0; j<=ymax_; j++) {
              box_[i][j] = new std::vector<int>[zmax_+1];
            }
        }
      // put the spheres in the boxes
      for (i=0; i<n_spheres; i++) {
          int ixmin, iymin, izmin, ixmax, iymax, izmax;
          SCVector3 l(sphere[i].center()), u(sphere[i].center());
          for (int j=0; j<3; j++) {
              l[j] -= probe_r + sphere[i].radius();
              u[j] += probe_r + sphere[i].radius();
            }
          get_box(l,ixmin,iymin,izmin);
          get_box(u,ixmax,iymax,izmax);
          for (int ii=ixmin; ii<=ixmax; ii++) {
              for (int jj=iymin; jj<=iymax; jj++) {
                  for (int kk=izmin; kk<=izmax; kk++) {
                      box_[ii][jj][kk].push_back(i);
                    }
                }
            }
        }
    }
}

int
ConnollyShape::get_box(const SCVector3 &v, int &x, int &y, int &z) const
{
  if (!box_) return 0;
  SCVector3 pos = v-lower_;
  x = (int)(pos[0]/l_);
  y = (int)(pos[1]/l_);
  z = (int)(pos[2]/l_);
  if (x<0) x=0;
  if (y<0) y=0;
  if (z<0) z=0;
  if (x>xmax_) x=xmax_;
  if (y>ymax_) y=ymax_;
  if (z>zmax_) z=zmax_;
  return 1;
}

ConnollyShape::~ConnollyShape()
{
  clear();
}

void
ConnollyShape::clear()
{
  delete[] sphere;
  sphere = 0;
  if (box_) {
      for (int i=0; i<=xmax_; i++) {
          for (int j=0; j<=ymax_; j++) {
              delete[] box_[i][j];
            }
          delete[] box_[i];
        }
      delete[] box_;
      box_ = 0;
    }
}

double
ConnollyShape::distance_to_surface(const SCVector3&r, SCVector3*grad) const
{
#if COUNT_CONNOLLY
  n_total_++;
#endif
  
  // can't compute grad so zero it if it is requested
  if (grad) {
      *grad = 0.0;
    }

  CS2Sphere probe_centers(r,probe_r);

  const int max_local_spheres = 60;
  CS2Sphere local_sphere[max_local_spheres];

  const double outside = 1.0;
  const double inside = -1.0;

  // find out which spheres are near the probe_centers sphere
  int n_local_spheres = 0;
  int boxi, boxj, boxk;
  if (get_box(r,boxi,boxj,boxk)) {
      std::vector<int> & box = box_[boxi][boxj][boxk];
      for (int ibox=0; ibox<box.size(); ibox++) {
          int i = box[ibox];
          double distance = sphere[i].distance(probe_centers);
          double r_i = sphere[i].radius();
          if (distance < r_i + probe_r) {
              if (distance < r_i - probe_r) {
#if COUNT_CONNOLLY
                  n_inside_vdw_++;
#endif
                  return inside;
                }
              if (n_local_spheres == max_local_spheres) {
                  ExEnv::err0() << indent
                       << "ConnollyShape::distance_to_surface:"
                       << " max_local_spheres exceeded\n";
                  abort();
                }
              local_sphere[n_local_spheres] = sphere[i];
              n_local_spheres++;
            }
        }
    }

#if COUNT_CONNOLLY
  if (n_local_spheres >= CONNOLLYSHAPE_N_WITH_NSPHERE_DIM) {
      n_with_nsphere_[CONNOLLYSHAPE_N_WITH_NSPHERE_DIM-1]++;
    }
  else {
      n_with_nsphere_[n_local_spheres]++;
    }
#endif

  if (probe_centers.intersect(local_sphere,n_local_spheres)
      == 1) return inside;
  return outside;
}

void
ConnollyShape::boundingbox(double valuemin,
                            double valuemax,
                            SCVector3& p1, SCVector3& p2)
{
  int i,j;
  if (valuemin < -1.0 || valuemax > 1.0) {
      ExEnv::err0() << indent
           << "ConnollyShape::boundingbox: value out of range\n";
      abort();
    }

  if (n_spheres == 0) {
      for (i=0; i<3; i++) {
          p1[i] = 0.0;
          p2[i] = 0.0;
        }