molecule.cc

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

      for (int j=i; j<i+5 && j<natom(); j++) {
          os << scprintf(" %10.5f", mass(j));
        }
      os << endl;
    }
}

int
Molecule::atom_label_to_index(const char *l) const
{
  int i;
  for (i=0; i<natom(); i++) {
      if (label(i) && !strcmp(l,label(i))) return i;
    }
  return -1;
}

double*
Molecule::charges() const
{
  double*result = new double[natoms_];
  if (charges_) {
      memcpy(result, charges_, sizeof(double)*natom());
    }
  else {
      for (int i=0; i<natom(); i++) result[i] = Z_[i];
    }
  return result;
}

double
Molecule::charge(int iatom) const
{
  if (charges_) return charges_[iatom];
  return Z_[iatom];
}

double
Molecule::nuclear_charge() const
{
  int i;
  double c = 0.0;
  for (i=0; i<natom(); i++) {
      c += charge(i);
    }
  return c;
}

void Molecule::save_data_state(StateOut& so)
{
  so.put(natoms_);
  SavableState::save_state(pg_.pointer(),so);
  SavableState::save_state(geometry_units_.pointer(),so);
  SavableState::save_state(atominfo_.pointer(),so);
  if (natoms_) {
      so.put(Z_, natoms_);
      so.put_array_double(r_[0], natoms_*3);
      so.put(charges_,natoms_);
    }
  if (mass_) {
      so.put(1);
      so.put_array_double(mass_, natoms_);
    }
  else {
      so.put(0);
    }
  if (labels_){
      so.put(1);
      for (int i=0; i<natoms_; i++) {
          so.putstring(labels_[i]);
        }
    }
  else {
      so.put(0);
    }
}

Molecule::Molecule(StateIn& si):
  SavableState(si),
  natoms_(0), r_(0), Z_(0), mass_(0), labels_(0)
{
  if (si.version(::class_desc<Molecule>()) < 4) {
      ExEnv::errn() << "Molecule: cannot restore from old molecules" << endl;
      abort();
    }
  si.get(natoms_);
  pg_ << SavableState::restore_state(si);
  geometry_units_ << SavableState::restore_state(si);
  atominfo_ << SavableState::restore_state(si);
  if (natoms_) {
      si.get(Z_);
      r_ = new double*[natoms_];
      r_[0] = new double[natoms_*3];
      si.get_array_double(r_[0],natoms_*3);
      for (int i=1; i<natoms_; i++) {
          r_[i] = &(r_[0][i*3]);
        }
      if (si.version(::class_desc<Molecule>()) > 4) {
          si.get(charges_);
        }
      else {
          charges_ = 0;
        }
    }
  int test;
  si.get(test);
  if (test) {
      mass_ = new double[natoms_];
      si.get_array_double(mass_, natoms_);
    }
  si.get(test);
  if (test){
      labels_ = new char*[natoms_];
      for (int i=0; i<natoms_; i++) {
          si.getstring(labels_[i]);
        }
    }

  nuniq_ = 0;
  equiv_ = 0;
  nequiv_ = 0;
  atom_to_uniq_ = 0;
  init_symmetry_info();
}

int
Molecule::atom_to_unique_offset(int iatom) const
{
  int iuniq = atom_to_uniq_[iatom];
  int nequiv = nequiv_[iuniq];
  for (int i=0; i<nequiv; i++) {
      if (equiv_[iuniq][i] == iatom) return i;
    }
  ExEnv::errn() << "Molecule::atom_to_unique_offset: internal error"
               << endl;
  return -1;
}

void
Molecule::set_point_group(const Ref<PointGroup>&ppg, double tol)
{
  ExEnv::out0() << indent
       << "Molecule: setting point group to " << ppg->symbol()
       << endl;
  pg_ = new PointGroup(*ppg.pointer());

  double r[3];
  for (int i=0; i<3; i++) {
      r[i] = -pg_->origin()[i];
      pg_->origin()[i] = 0;
    }
  translate(r);

  clear_symmetry_info();
  init_symmetry_info();

  cleanup_molecule(tol);
}

Ref<PointGroup>
Molecule::point_group() const
{
  return pg_;
}

SCVector3
Molecule::center_of_mass() const
{
  SCVector3 ret;
  double M;

  ret = 0.0;
  M = 0.0;

  for (int i=0; i < natom(); i++) {
    double m = mass(i);
    ret += m * SCVector3(r(i));
    M += m;
  }

  ret *= 1.0/M;

  return ret;
}

double
Molecule::nuclear_repulsion_energy()
{
  int i, j;
  double e=0.0;

  for (i=1; i < natoms_; i++) {
    SCVector3 ai(r(i));
    double Zi = charge(i);
    
    for (j=0; j < i; j++) {
        SCVector3 aj(r(j));
        e += Zi * charge(j) / ai.dist(aj);
      }
    }

  return e;
}

void
Molecule::nuclear_repulsion_1der(int center, double xyz[3])
{
  int i,j,k;
  double rd[3],r2;
  double factor;

  xyz[0] = 0.0;
  xyz[1] = 0.0;
  xyz[2] = 0.0;
  for (i=0; i < natoms_; i++) {
      SCVector3 ai(r(i));
      double Zi = charge(i);

      for (j=0; j < i; j++) {
          if (center==i || center==j) {
              SCVector3 aj(r(j));

              r2 = 0.0;
              for (k=0; k < 3; k++) {
                  rd[k] = ai[k] - aj[k];
                  r2 += rd[k]*rd[k];
                }
        
              factor = - Zi * charge(j) * pow(r2,-1.5);
              if (center==j) factor = -factor;
              for (k=0; k<3; k++) {
                  xyz[k] += factor * rd[k];
                }
            }
        }
    }
}

void
Molecule::nuclear_charge_efield(const double *charges,
                                const double *position, double *efield)
{
  int i,j;
  double tmp;
  double rd[3];

  for (i=0; i<3; i++) efield[i] = 0.0;

  for (i=0; i<natoms_; i++) {
      SCVector3 a(r(i));
      tmp = 0.0;
      for (j=0; j<3; j++) {
          rd[j] = position[j] - a[j];
          tmp += rd[j]*rd[j];
        }
      tmp = charges[i]/(tmp*sqrt(tmp));
      for (j=0; j<3; j++) {
          efield[j] +=  rd[j] * tmp;
        }
    }
}

void
Molecule::nuclear_efield(const double *position, double *efield)
{
  int i,j;
  double tmp;
  double rd[3];

  for (i=0; i<3; i++) efield[i] = 0.0;

  for (i=0; i<natoms_; i++) {
      SCVector3 a(r(i));
      tmp = 0.0;
      for (j=0; j<3; j++) {
          rd[j] = position[j] - a[j];
          tmp += rd[j]*rd[j];
        }
      tmp = charge(i)/(tmp*sqrt(tmp));
      for (j=0; j<3; j++) {
          efield[j] +=  rd[j] * tmp;
        }
    }
}

int
Molecule::atom_at_position(double *v, double tol) const
{
  SCVector3 p(v);
  for (int i=0; i < natom(); i++) {
      SCVector3 ai(r(i));
      if (p.dist(ai) < tol) return i;
    }
  return -1;
}

void
Molecule::symmetrize(const Ref<PointGroup> &pg, double tol)
{
  pg_ = new PointGroup(pg);

  // translate to the origin of the symmetry frame
  double r[3];
  for (int i=0; i<3; i++) {
      r[i] = -pg_->origin()[i];
      pg_->origin()[i] = 0;
    }
  translate(r);

  symmetrize(tol);
}

// We are given a molecule which may or may not have just the symmetry
// distinct atoms in it.  We have to go through the existing set of atoms,
// perform each symmetry operation in the point group on each of them, and
// then add the new atom if it isn't in the list already

void
Molecule::symmetrize(double tol)
{
  // if molecule is c1, don't do anything
  if (!strcmp(this->point_group()->symbol(),"c1")) {
    init_symmetry_info();
    return;
    }

  clear_symmetry_info();

  Molecule *newmol = new Molecule(*this);
  
  CharacterTable ct = this->point_group()->char_table();

  SCVector3 np;
  SymmetryOperation so;
  
  for (int i=0; i < natom(); i++) {
    SCVector3 ac(r(i));

    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];
      }

      int atom = newmol->atom_at_position(np.data(), tol);
      if (atom < 0) {
        newmol->add_atom(Z_[i],np[0],np[1],np[2],label(i));
      }
      else {
        if (Z(i) != newmol->Z(atom)
            || fabs(mass(i)-newmol->mass(atom))>1.0e-10) {
            ExEnv::err0() << "Molecule: symmetrize: atom mismatch" << endl;
            abort();
        }
      }
    }
  }
  
  Ref<Units> saved_units = geometry_units_;
  *this = *newmol;
  geometry_units_ = saved_units;
  delete newmol;

  init_symmetry_info();
}

void
Molecule::translate(const double *r)
{
  for (int i=0; i < natom(); i++) {
    r_[i][0] += r[0];
    r_[i][1] += r[1];
    r_[i][2] += r[2];
  }
}

// move the molecule to the center of mass
void
Molecule::move_to_com()
{
  SCVector3 com = -center_of_mass();
  translate(com.data());
}

// find the 3 principal coordinate axes, and rotate the molecule to be 
// aligned along them.  also rotate the symmetry frame contained in point_group
void
Molecule::transform_to_principal_axes(int trans_frame)
{
  // mol_move_to_com(mol);

  double *inert[3], inert_dat[9], *evecs[3], evecs_dat[9];
  double evals[3];

  int i,j,k;
  for (i=0; i < 3; i++) {
    inert[i] = &inert_dat[i*3];
    evecs[i] = &evecs_dat[i*3];
  }
  memset(inert_dat,'\0',sizeof(double)*9);
  memset(evecs_dat,'\0',sizeof(double)*9);

  for (i=0; i < natom(); i++) {
    SCVector3 ac(r(i));
    double m=mass(i);
    inert[0][0] += m * (ac[1]*ac[1] + ac[2]*ac[2]);
    inert[1][0] -= m * ac[0]*ac[1];
    inert[1][1] += m * (ac[0]*ac[0] + ac[2]*ac[2]);
    inert[2][0] -= m * ac[0]*ac[2];
    inert[2][1] -= m * ac[1]*ac[2];
    inert[2][2] += m * (ac[0]*ac[0] + ac[1]*ac[1]);
  }

  inert[0][1] = inert[1][0];
  inert[0][2] = inert[2][0];
  inert[1][2] = inert[2][1];

 // cleanup inert
  for (i=0; i < 3; i++) {
    for (int j=0; j <= i; j++) {
      if (fabs(inert[i][j]) < 1.0e-5) {
        inert[i][j]=inert[j][i]=0.0;
      }