molsymm.cc

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

//
// molsymm.cc
//
// Copyright (C) 1997 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.
//

#include <math.h>

#include <util/misc/formio.h>
#include <chemistry/molecule/molecule.h>

using namespace std;
using namespace sc;

#undef DEBUG

void
Molecule::clear_symmetry_info()
{
  for (int i=0; i<nuniq_; i++) {
    delete[] equiv_[i];
    }
  delete[] equiv_;
  delete[] nequiv_;
  delete[] atom_to_uniq_;
  nuniq_ = 0;
  equiv_ = 0;
  nequiv_ = 0;
  atom_to_uniq_ = 0;
}

void
Molecule::init_symmetry_info(double tol)
{
  if (equiv_)
    clear_symmetry_info();
  
  if (natom() == 0) {
    nuniq_ = 0;
    equiv_ = 0;
    nequiv_ = 0;
    atom_to_uniq_ = 0;
    return;
    }

  nequiv_ = new int[natom()];
  atom_to_uniq_ = new int[natom()];
  equiv_ = new int*[natom()];

  if (!strcmp(point_group()->symbol(),"c1")) {
    nuniq_ = natom();
    for (int i=0; i < natom(); i++) {
      nequiv_[i]=1;
      equiv_[i]=new int[1];
      equiv_[i][0]=i;
      atom_to_uniq_[i]=i;
      }
    return;
  }

  // the first atom is always unique
  nuniq_ = 1;
  nequiv_[0]=1;
  equiv_[0] = new int[1];
  equiv_[0][0]=0;
  atom_to_uniq_[0]=0;

  CharacterTable ct = point_group()->char_table();

  SCVector3 ac;
  SymmetryOperation so;
  SCVector3 np;

  // find the equivalent atoms
  int i;
  for (i=1; i < natom(); i++) {
    ac = r(i);
    int i_is_unique=1;
    int i_equiv=0;

    // apply all symmetry ops in the group to the atom
    for (int g=0; g < ct.order(); g++) {
      so = ct.symm_operation(g);
      for (int ii=0; ii < 3; ii++) {
        np[ii]=0;
        for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * ac[jj];
        }

      // see if the transformed atom is equivalent to a unique atom
      for (int j=0; j<nuniq_; j++) {
        int uniq = equiv_[j][0];
        SCVector3 aj(r(uniq));
        if (np.dist(aj) < tol
            && Z(uniq) == Z(i)
            && fabs(charge(uniq)-charge(i)) < tol
            && fabs(mass(uniq)-mass(i)) < tol) {
          i_is_unique = 0;
          i_equiv = j;
          break;
          }
        }
      }
    if (i_is_unique) {
      nequiv_[nuniq_]=1;
      equiv_[nuniq_]=new int[1];
      equiv_[nuniq_][0]=i;
      atom_to_uniq_[i] = nuniq_;
      nuniq_++;
      }
    else {
      int *tmp = new int[nequiv_[i_equiv]+1];
      memcpy(tmp,equiv_[i_equiv],nequiv_[i_equiv]*sizeof(int));
      delete[] equiv_[i_equiv];
      equiv_[i_equiv] = tmp;
      equiv_[i_equiv][nequiv_[i_equiv]] = i;
      nequiv_[i_equiv]++;
      atom_to_uniq_[i] = i_equiv;
      }
    }

  // The first atom in the equiv list is considered the primary unique
  // atom.  Just to make things look pretty, make the atom with the most
  // zeros in its x, y, z coordinate the unique atom.  Nothing else should
  // rely on this being done.
  double ztol=1.0e-5;
  for (i=0; i < nuniq_; i++) {
    int maxzero = 0;
    int jmaxzero = 0;
    for (int j=0; j<nequiv_[i]; j++) {
      int nzero = 0;
      for (int k=0; k<3; k++) if (fabs(r(equiv_[i][j],k)) < ztol) nzero++;
      if (nzero > maxzero) {
        maxzero = nzero;
        jmaxzero = j;
        }
      }
    int tmp = equiv_[i][jmaxzero];
    equiv_[i][jmaxzero] = equiv_[i][0];
    equiv_[i][0] = tmp;
    }
}

int
Molecule::has_inversion(SCVector3 &origin, double tol) const
{
  for (int i=0; i<natom(); i++) {
    SCVector3 inverted = origin-(SCVector3(r(i))-origin);
    int atom = atom_at_position(inverted.data(), tol);
    if (atom < 0
        || Z(atom) != Z(i)
        || fabs(charge(atom)-charge(i)) > tol
        || fabs(mass(atom)-mass(i)) > tol) {
      return 0;
      }
    }
  return 1;
}

int
Molecule::is_plane(SCVector3 &origin, SCVector3 &uperp, double tol) const
{
  for (int i=0; i<natom(); i++) {
    SCVector3 A = SCVector3(r(i))-origin;
    SCVector3 Apar = uperp.dot(A) * uperp;
    SCVector3 Aperp = A - Apar;
    A = (Aperp - Apar) + origin;
    int atom = atom_at_position(A.data(), tol);
    if (atom < 0
        || Z(atom) != Z(i)
        || fabs(charge(atom)-charge(i)) > tol
        || fabs(mass(atom)-mass(i)) > tol) {
      //ExEnv::outn() << "  is_plane: rejected (atom " << i << ")" << endl;
      return 0;
      }
    }
  return 1;
}

int
Molecule::is_axis(SCVector3 &origin, SCVector3 &axis,
                  int order, double tol) const
{
  // loop through atoms to see if axis is a c2 axis
  for (int i=0; i<natom(); i++) {
    SCVector3 A = SCVector3(r(i))-origin;
    for (int j=1; j<order; j++) {
      SCVector3 R = A;
      R.rotate(j*2.0*M_PI/order, axis);
      R += origin;
      int atom = atom_at_position(R.data(), tol);
      if (atom < 0
          || Z(atom) != Z(i)
          || fabs(charge(atom)-charge(i)) > tol
          || fabs(mass(atom)-mass(i)) > tol) {
        //ExEnv::outn() << "  is_axis: rejected (atom " << i << ")" << endl;
        return 0;
        }
      }
    }
  return 1;
}

enum AxisName { XAxis, YAxis, ZAxis };

static AxisName
like_world_axis(SCVector3 &axis,
                const SCVector3 &worldxaxis,
                const SCVector3 &worldyaxis,
                const SCVector3 &worldzaxis
                )
{
  AxisName like;
  double xlikeness = fabs(axis.dot(worldxaxis));
  double ylikeness = fabs(axis.dot(worldyaxis));
  double zlikeness = fabs(axis.dot(worldzaxis));
  if (xlikeness > ylikeness && xlikeness > zlikeness) {
    like = XAxis;
    if (axis.dot(worldxaxis) < 0) axis = - axis; 
    }
  else if (ylikeness > zlikeness) {
    like = YAxis;
    if (axis.dot(worldyaxis) < 0) axis = - axis; 
    }
  else {
    like = ZAxis;
    if (axis.dot(worldzaxis) < 0) axis = - axis; 
    }
  return like;
}

Ref<PointGroup>
Molecule::highest_point_group(double tol) const
{
  int i,j;

  SCVector3 com = center_of_mass();

  SCVector3 worldzaxis(0.,0.,1.);
  SCVector3 worldxaxis(1.,0.,0.);
  SCVector3 worldyaxis(0.,1.,0.);

  int linear,planar;
  is_linear_planar(linear,planar,tol);

  int have_inversion = has_inversion(com,tol);

  // check for C2 axis
  SCVector3 c2axis;
  int have_c2axis = 0;
  if (natom() < 2) {
    have_c2axis = 1;
    c2axis = SCVector3(0.0,0.0,1.0);
    }
  else if (linear) {
    have_c2axis = 1;
    c2axis = SCVector3(r(1)) - SCVector3(r(0));
    c2axis.normalize();
    }
  else if (planar && have_inversion) {
    // there is a c2 axis that won't be found using the usual algorithm.
    // find two noncolinear atom-atom vectors (we know that linear==0)
    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_c2axis = 1;
        BAxCA.normalize();
        c2axis = BAxCA;
        break;
        }
      }
    }
  else {
    // loop through pairs of atoms to find C2 axis candidates
    for (i=0; i<natom(); i++) {
      SCVector3 A = SCVector3(r(i))-com;
      double AdotA = A.dot(A);
      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 axis = A+B;
        // atoms colinear with the com don't work
        if (axis.norm() < tol) continue;
        axis.normalize();
        if (is_axis(com,axis,2,tol)) {
          have_c2axis = 1;
          c2axis = axis;
          goto found_c2axis;
          }
        }
      }
    }
  found_c2axis:

  AxisName c2like = ZAxis;
  if (have_c2axis) {
    // try to make the sign of the axis correspond to one of the world axes
    c2like = like_world_axis(c2axis,worldxaxis,worldyaxis,worldzaxis);
    }

  // check for C2 axis perp to first C2 axis
  SCVector3 c2axisperp;
  int have_c2axisperp = 0;
  if (have_c2axis) {
    if (natom() < 2) {
      have_c2axisperp = 1;
      c2axisperp = SCVector3(1.0,0.0,0.0);
      }
    else if (linear) {
      if (have_inversion) {
        have_c2axisperp = 1;
        c2axisperp = c2axis.perp_unit(SCVector3(0.0,0.0,1.0));
        }
      }
    else {
      // loop through pairs of atoms to find C2 axis candidates
      for (i=0; i<natom(); i++) {
        SCVector3 A = SCVector3(r(i))-com;
        double AdotA = A.dot(A);
        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 axis = A+B;
          // atoms colinear with the com don't work
          if (axis.norm() < tol) continue;
          axis.normalize();