molsymm.cc

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

          // if axis is not perp continue
          if (fabs(axis.dot(c2axis)) > tol) continue;
          if (is_axis(com,axis,2,tol)) {
            have_c2axisperp = 1;
            c2axisperp = axis;
            goto found_c2axisperp;
            }
          }
        }
      }
    }
  found_c2axisperp:

  AxisName c2perplike;
  if (have_c2axisperp) {
    // try to make the sign of the axis correspond to one of the world axes
    c2perplike = like_world_axis(c2axisperp,worldxaxis,worldyaxis,worldzaxis);

    // try to make c2axis the z axis
    if (c2perplike == ZAxis) {
      SCVector3 tmpv = c2axisperp;
      tmpv = c2axisperp; c2axisperp = c2axis; c2axis = tmpv;
      c2perplike = c2like;
      c2like = ZAxis;
      }
    if (c2like != ZAxis) {
      if (c2like == XAxis) c2axis = c2axis.cross(c2axisperp);
      else c2axis = c2axisperp.cross(c2axis);
      c2like = like_world_axis(c2axis,worldxaxis,worldyaxis,worldzaxis);
      }

    // try to make c2axisperp like the x axis
    if (c2perplike == YAxis) {
      c2axisperp = c2axisperp.cross(c2axis);
      c2perplike = like_world_axis(c2axisperp,
                                   worldxaxis,worldyaxis,worldzaxis);
      }
    }

  // check for vertical plane
  int have_sigmav = 0;
  SCVector3 sigmav;
  if (have_c2axis) {
    if (natom() < 2) {
      have_sigmav = 1;
      sigmav = c2axisperp;
      }
    else if (linear) {
      have_sigmav = 1;
      if (have_c2axisperp) {
        sigmav = c2axisperp;
        }
      else {
        sigmav = c2axis.perp_unit(SCVector3(0.0,0.0,1.0));
        }
      }
    else {
      // loop through pairs of atoms to find sigma v plane candidates
      for (i=0; i<natom(); i++) {
        SCVector3 A = SCVector3(r(i))-com;
        double AdotA = A.dot(A);
        // the second atom can equal i because i might be in the plane.
        for (j=0; j<=i; j++) {
          // the atoms must be identical
          if (Z(i) != Z(j) || fabs(mass(i)-mass(j)) > tol) continue;
          SCVector3 B = SCVector3(r(j))-com;
          // the atoms must be the same distance from the com
          if (fabs(AdotA - B.dot(B)) > tol) continue;
          SCVector3 inplane = B+A;
          double norm_inplane = inplane.norm();
          if (norm_inplane < tol) continue;
          inplane *= 1.0/norm_inplane;
          SCVector3 perp = c2axis.cross(inplane);
          double norm_perp = perp.norm();
          if (norm_perp < tol) continue;
          perp *= 1.0/norm_perp;
          if (is_plane(com,perp,tol)) {
            have_sigmav = 1;
            sigmav = perp;
            goto found_sigmav;
            }
          }
        }
      }
    }
  found_sigmav:

  if (have_sigmav) {
    // try to make the sign of the oop vec correspond to one of the world axes
    int sigmavlike = like_world_axis(sigmav,worldxaxis,worldyaxis,worldzaxis);

    // choose sigmav to be the world x axis, if possible
    if (c2like == ZAxis && sigmavlike == YAxis) {
      sigmav = sigmav.cross(c2axis);
      }
    else if (c2like == YAxis && sigmavlike == ZAxis) {
      sigmav = c2axis.cross(sigmav);
      }
    }

  // under certain conditions i need to know if there is any sigma plane
  int have_sigma = 0;
  SCVector3 sigma;
  if (!have_inversion && !have_c2axis) {
    if (planar) {
      // find two noncolinear atom-atom vectors
      // we know that linear==0 since !have_c2axis
      SCVector3 BA = SCVector3(r(1))-SCVector3(r(0));
      BA.normalize();
      for (i=2; i<natom(); i++) {
        SCVector3 CA = SCVector3(r(i))-SCVector3(r(0));
        CA.normalize();
        SCVector3 BAxCA = BA.cross(CA);
        if (BAxCA.norm() > tol) {
          have_sigma = 1;
          BAxCA.normalize();
          sigma = BAxCA;
          break;
          }
        }
      }
    else {
      // loop through pairs of atoms to construct trial planes
      for (i=0; i<natom(); i++) {
        SCVector3 A = SCVector3(r(i))-com;
        double AdotA = A.dot(A);
        for (j=0; j<i; j++) {
          //ExEnv::outn() << "sigma atoms = " << i << ", " << j << endl;
          // the atoms must be identical
          if (Z(i) != Z(j) || fabs(mass(i)-mass(j)) > tol) continue;
          SCVector3 B = SCVector3(r(j))-com;
          double BdotB = B.dot(B);
          // the atoms must be the same distance from the com
          if (fabs(AdotA - BdotB) > tol) continue;
          SCVector3 perp = B-A;
          double norm_perp = perp.norm();
          if (norm_perp < tol) {
            //ExEnv::outn() << "  rejected (atoms at same point?)" << endl;
            continue;
            }
          perp *= 1.0/norm_perp;
          if (is_plane(com,perp,tol)) {
            have_sigma = 1;
            sigma = perp;
            goto found_sigma;
            }
          }
        }
      }
    }
  found_sigma:

  if (have_sigma) {
    // try to make the sign of the oop vec correspond to one of the world axes
    double xlikeness = fabs(sigma.dot(worldxaxis));
    double ylikeness = fabs(sigma.dot(worldyaxis));
    double zlikeness = fabs(sigma.dot(worldzaxis));
    if (xlikeness > ylikeness && xlikeness > zlikeness) {
      if (sigma.dot(worldxaxis) < 0) sigma = - sigma; 
      }
    else if (ylikeness > zlikeness) {
      if (sigma.dot(worldyaxis) < 0) sigma = - sigma; 
      }
    else {
      if (sigma.dot(worldzaxis) < 0) sigma = - sigma; 
      }
    }

#ifdef DEBUG
  ExEnv::out0()
       << indent << "highest point group:" << endl
       << indent << "  linear          = " << linear << endl
       << indent << "  planar          = " << planar << endl
       << indent << "  have_inversion  = " << have_inversion << endl
       << indent << "  have_c2axis     = " << have_c2axis << endl
       << indent << "  have_c2axisperp = " << have_c2axisperp << endl
       << indent << "  have_sigmav     = " << have_sigmav << endl
       << indent << "  have_sigma      = " << have_sigma << endl;

  if (have_c2axis)
    ExEnv::out0() << indent << "  c2axis      = " << c2axis << endl;
  if (have_c2axisperp)
    ExEnv::out0() << indent << "  c2axisperp  = " << c2axisperp << endl;
  if (have_sigmav)
    ExEnv::out0() << indent << "  sigmav      = " << sigmav << endl;
  if (have_sigma)
    ExEnv::out0() << indent << "  sigma       = " << sigma << endl;
#endif

  // Find the three axes for the symmetry frame
  SCVector3 xaxis = worldxaxis;
  SCVector3 yaxis;
  SCVector3 zaxis = worldzaxis;;
  if (have_c2axis) {
    zaxis = c2axis;
    if (have_sigmav) {
      xaxis = sigmav;
      }
    else if (have_c2axisperp) {
      xaxis = c2axisperp;
      }
    else {
      // any axis orthogonal to the zaxis will do
      xaxis = zaxis.perp_unit(zaxis);
      }
    }
  else if (have_sigma) {
    zaxis = sigma;
    xaxis = zaxis.perp_unit(zaxis);
    }
  // the y axis is then -x cross z
  yaxis = - xaxis.cross(zaxis);

#ifdef DEBUG
  ExEnv::outn() << "X: " << xaxis << endl;
  ExEnv::outn() << "Y: " << yaxis << endl;
  ExEnv::outn() << "Z: " << zaxis << endl;
#endif

  SymmetryOperation frame;
  SCVector3 origin;
  for (i=0; i<3; i++) {
    frame(i,0) = xaxis[i];
    frame(i,1) = yaxis[i];
    frame(i,2) = zaxis[i];
    origin[i] = com[i];
    }

#ifdef DEBUG
  ExEnv::out0() << "frame:" << endl;
  frame.print(ExEnv::out0());

  ExEnv::out0() << "origin:" << endl;
  origin.print(ExEnv::out0());
#endif

  Ref<PointGroup> pg;
  if (have_inversion) {
    if (have_c2axis) {
      if (have_sigmav) {
        pg = new PointGroup("d2h",frame,origin);
        }
      else {
        pg = new PointGroup("c2h",frame,origin);
        }
      }
    else {
      pg = new PointGroup("ci",frame,origin);
      }
    }
  else {
    if (have_c2axis) {
      if (have_sigmav) {
        pg = new PointGroup("c2v",frame,origin);
        }
      else {
        if (have_c2axisperp) {
          pg = new PointGroup("d2",frame,origin);
          }
        else {
          pg = new PointGroup("c2",frame,origin);
          }
        }
      }
    else {
      if (have_sigma) {
        pg = new PointGroup("cs",frame,origin);
        }
      else {
        pg = new PointGroup("c1",frame,origin);
        }
      }
    }

  return pg;
}

int
Molecule::is_linear(double tol) const
{
  int linear, planar;

  is_linear_planar(linear,planar,tol);

  return linear;
}

int
Molecule::is_planar(double tol) const
{
  int linear, planar;

  is_linear_planar(linear,planar,tol);

  return planar;
}

void
Molecule::is_linear_planar(int &linear, int &planar, double tol) const
{
  if (natom() < 3) {
    linear = 1;
    planar = 1;
    return;
    }

  // find three atoms not on the same line
  SCVector3 A = r(0);
  SCVector3 B = r(1);
  SCVector3 BA = B-A;
  BA.normalize();
  SCVector3 CA;

  int i;
  double min_BAdotCA = 1.0;
  for (i=2; i<natom(); i++) {
    SCVector3 tmp = SCVector3(r(i))-A;
    tmp.normalize();
    if (fabs(BA.dot(tmp)) < min_BAdotCA) {
      CA = tmp;
      min_BAdotCA = fabs(BA.dot(tmp));
      }
    }
  if (min_BAdotCA >= 1.0 - tol) {
    linear = 1;
    planar = 1;
    return;
    }

  linear = 0;
  if (natom() < 4) {
    planar = 1;
    return;
    }

  // check for nontrivial planar molecules
  SCVector3 BAxCA = BA.cross(CA);
  BAxCA.normalize();
  for (i=2; i<natom(); i++) {
    SCVector3 tmp = SCVector3(r(i))-A;
    if (fabs(tmp.dot(BAxCA)) > tol) {
      planar = 0;
      return;
      }
    }
  planar = 1;
  return;
}

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

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