molecule.cc

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

//
// molecule.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 <math.h>
#include <string.h>
#include <stdio.h>

#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <chemistry/molecule/molecule.h>
#include <chemistry/molecule/formula.h>
#include <chemistry/molecule/localdef.h>
#include <math/scmat/cmatrix.h>

using namespace std;
using namespace sc;

//////////////////////////////////////////////////////////////////////////
// Molecule

static ClassDesc Molecule_cd(
  typeid(Molecule),"Molecule",5,"public SavableState",
  create<Molecule>, create<Molecule>, create<Molecule>);

Molecule::Molecule():
  natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
{
  pg_ = new PointGroup;
  atominfo_ = new AtomInfo();
  geometry_units_ = new Units("bohr");
  nuniq_ = 0;
  equiv_ = 0;
  nequiv_ = 0;
  atom_to_uniq_ = 0;
  init_symmetry_info();
}

Molecule::Molecule(const Molecule& mol):
 natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
{
  nuniq_ = 0;
  equiv_ = 0;
  nequiv_ = 0;
  atom_to_uniq_ = 0;
  *this=mol;
}

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

void
Molecule::clear()
{
  if (r_) {
      delete[] r_[0];
      delete[] r_;
      r_ = 0;
    }
  if (labels_) {
      for (int i=0; i<natoms_; i++) {
          delete[] labels_[i];
        }
      delete[] labels_;
      labels_ = 0;
    }
  delete[] charges_;
  charges_ = 0;
  delete[] mass_;
  mass_ = 0;
  delete[] Z_;
  Z_ = 0;

  clear_symmetry_info();
}

Molecule::Molecule(const Ref<KeyVal>&input):
 natoms_(0), r_(0), Z_(0), charges_(0), mass_(0), labels_(0)
{
  nuniq_ = 0;
  equiv_ = 0;
  nequiv_ = 0;
  atom_to_uniq_ = 0;

  atominfo_ << input->describedclassvalue("atominfo");
  if (atominfo_.null()) atominfo_ = new AtomInfo;
  if (input->exists("pdb_file")) {
      geometry_units_ = new Units("angstrom");
      char* filename = input->pcharvalue("pdb_file");
      read_pdb(filename);
      delete[] filename;
    }
  else {
      // check for old style units input first
      if (input->booleanvalue("angstrom")
          ||input->booleanvalue("angstroms")
          ||input->booleanvalue("aangstrom")
          ||input->booleanvalue("aangstroms")) {
          geometry_units_ = new Units("angstrom");
        }
      // check for new style units input
      else {
          geometry_units_ = new Units(input->pcharvalue("unit"),
                                      Units::Steal);
        }
        
      double conv = geometry_units_->to_atomic_units();

      // get the number of atoms and make sure that the geometry and the
      // atoms array have the same number of atoms.
      // right now we read in the unique atoms...then we will symmetrize.
      // the length of atoms must still equal the length of geometry, but
      // we'll try to set up atom_labels such that different lengths are
      // possible
      int natom = input->count("geometry");
      if (natom != input->count("atoms")) {
          ExEnv::err0() << indent
               << "Molecule: size of \"geometry\" != size of \"atoms\"\n";
          abort();
        }

      int i;
      for (i=0; i<natom; i++) {
          char *name, *label;
          int ghost = input->booleanvalue("ghost",i);
          double charge = input->doublevalue("charge",i);
          int have_charge = input->error() == KeyVal::OK;
          if (ghost) {
              have_charge = 1;
              charge = 0.0;
            }
          add_atom(AtomInfo::string_to_Z(name = input->pcharvalue("atoms",i)),
                   input->doublevalue("geometry",i,0)*conv,
                   input->doublevalue("geometry",i,1)*conv,
                   input->doublevalue("geometry",i,2)*conv,
                   label = input->pcharvalue("atom_labels",i),
                   input->doublevalue("mass",i),
                   have_charge, charge
              );
          delete[] name;
          delete[] label;
        }
    }

  char *symmetry = input->pcharvalue("symmetry");
  double symtol = input->doublevalue("symmetry_tolerance",
                                     KeyValValuedouble(1.0e-4));
  nuniq_ = 0;
  equiv_ = 0;
  nequiv_ = 0;
  atom_to_uniq_ = 0;
  if (symmetry && !strcmp(symmetry,"auto")) {
      set_point_group(highest_point_group(symtol), symtol*10.0);
    }
  else {
      pg_ = new PointGroup(input);

      // 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);

      if (input->booleanvalue("redundant_atoms")) {
          init_symmetry_info();
          cleanup_molecule(symtol);
        }
      else {
          symmetrize();
          // In case we were given redundant atoms, clean up
          // the geometry so the symmetry is exact.
          cleanup_molecule(symtol);
        }
    }
  delete[] symmetry;
}

Molecule&
Molecule::operator=(const Molecule& mol)
{
  clear();

  pg_ = new PointGroup(*(mol.pg_.pointer()));
  atominfo_ = mol.atominfo_;
  geometry_units_ = new Units(mol.geometry_units_->string_rep());

  natoms_ = mol.natoms_;

  if (natoms_) {
      if (mol.mass_) {
          mass_ = new double[natoms_];
          memcpy(mass_,mol.mass_,natoms_*sizeof(double));
        }
      if (mol.charges_) {
          charges_ = new double[natoms_];
          memcpy(charges_,mol.charges_,natoms_*sizeof(double));
        }
      if (mol.labels_) {
          labels_ = new char *[natoms_];
          for (int i=0; i<natoms_; i++) {
              if (mol.labels_[i]) {
                  labels_[i] = strcpy(new char[strlen(mol.labels_[i])+1],
                                      mol.labels_[i]);
                }
              else labels_[i] = 0;
            }
        }
      r_ = new double*[natoms_];
      r_[0] = new double[natoms_*3];
      for (int i=0; i<natoms_; i++) {
          r_[i] = &(r_[0][i*3]);
        }
      memcpy(r_[0], mol.r_[0], natoms_*3*sizeof(double));
      Z_ = new int[natoms_];
      memcpy(Z_, mol.Z_, natoms_*sizeof(int));
    }

  init_symmetry_info();

  return *this;
}

void
Molecule::add_atom(int Z,double x,double y,double z,
                   const char *label,double mass,
                   int have_charge, double charge)
{
  int i;

  // allocate new arrays
  int *newZ = new int[natoms_+1];
  double **newr = new double*[natoms_+1];
  double *newr0 = new double[(natoms_+1)*3];
  char **newlabels = 0;
  if (label || labels_) {
      newlabels = new char*[natoms_+1];
    }
  double *newcharges = 0;
  if (have_charge || charges_) {
      newcharges = new double[natoms_+1];
    }
  double *newmass = 0;
  if (mass_ || mass != 0.0) {
      newmass = new double[natoms_+1];
    }

  // setup the r_ pointers
  for (i=0; i<=natoms_; i++) {
      newr[i] = &(newr0[i*3]);
    }

  // copy old data to new arrays
  if (natoms_) {
      memcpy(newZ,Z_,sizeof(int)*natoms_);
      memcpy(newr0,r_[0],sizeof(double)*natoms_*3);
      if (labels_) {
          memcpy(newlabels,labels_,sizeof(char*)*natoms_);
        }
      else if (newlabels) {
          memset(newlabels,0,sizeof(char*)*natoms_);
        }
      if (charges_) {
          memcpy(newcharges,charges_,sizeof(double)*natoms_);
        }
      else if (newcharges) {
          for (i=0; i<natoms_; i++) newcharges[i] = Z_[i];
        }
      if (mass_) {
          memcpy(newmass,mass_,sizeof(double)*natoms_);
        }
      else if (newmass) {
          memset(newmass,0,sizeof(double)*natoms_);
        }
    }

  // delete old data
  delete[] Z_;
  if (r_) {
      delete[] r_[0];
      delete[] r_;
    }
  delete[] labels_;
  delete[] charges_;
  delete[] mass_;

  // setup new pointers
  Z_ = newZ;
  r_ = newr;
  labels_ = newlabels;
  charges_ = newcharges;
  mass_ = newmass;

  // copy info for this atom into arrays
  Z_[natoms_] = Z;
  r_[natoms_][0] = x;
  r_[natoms_][1] = y;
  r_[natoms_][2] = z;
  if (mass_) mass_[natoms_] = mass;
  if (label) {
      labels_[natoms_] = strcpy(new char[strlen(label)+1],label);
    }
  else if (labels_) {
      labels_[natoms_] = 0;
    }
  if (have_charge) {
      charges_[natoms_] = charge;
    }
  else if (charges_) {
      charges_[natoms_] = Z;
    }

  natoms_++;
}

void
Molecule::print_parsedkeyval(ostream& os,
                             int print_pg,
                             int print_unit,
                             int number_atoms) const
{
  int i;

  double conv = geometry_units_->from_atomic_units();

  if (print_pg) pg_->print(os);
  if (print_unit && geometry_units_->string_rep()) {
      os << indent
         << "unit = \"" << geometry_units_->string_rep() << "\""
         << endl;
    }
  os << indent << "{";
  if (number_atoms) os << scprintf("%3s", "n");
  os << scprintf(" %5s", "atoms");
  if (labels_) os << scprintf(" %11s", "atom_labels");
  int int_charges = 1;
  if (charges_) {
      for (i=0;i<natom();i++) if (charges_[i]!=(int)charges_[i]) int_charges=0;
      if (int_charges) {
          os << scprintf(" %7s", "charge");
        }
      else {
          os << scprintf(" %17s", "charge");
        }
    }
  os << scprintf("  %16s", "")
     << scprintf(" %16s", "geometry   ")
     << scprintf(" %16s ", "");
  os << "}={" << endl;
  for (i=0; i<natom(); i++) {
      os << indent;
      if (number_atoms) os << scprintf(" %3d", i+1);
      os << scprintf(" %5s", AtomInfo::symbol(Z_[i]));
      if (labels_) {
          const char *lab = labels_[i];
          if (lab == 0) lab = "";
          char  *qlab = new char[strlen(lab)+3];
          strcpy(qlab,"\"");
          strcat(qlab,lab);
          strcat(qlab,"\"");
          os << scprintf(" %11s",qlab);
          delete[] qlab;
        }
      if (charges_) {
          if (int_charges) os << scprintf(" %7.4f", charges_[i]);
          else os << scprintf(" %17.15f", charges_[i]);
        }
      os << scprintf(" [% 16.10f", conv * r(i,0))
         << scprintf(" % 16.10f", conv * r(i,1))
         << scprintf(" % 16.10f]", conv * r(i,2))
         << endl;
    }
  os << indent << "}" << endl;
}

void
Molecule::print(ostream& os) const
{
  int i;

  MolecularFormula *mf = new MolecularFormula(this);
  os << indent
     << "Molecular formula: " << mf->formula() << endl;
  delete mf;

  os << indent << "molecule<Molecule>: (" << endl;
  os << incindent;
  print_parsedkeyval(os);
  os << decindent;
  os << indent << ")" << endl;

  os << indent << "Atomic Masses:" << endl;
  for (i=0; i<natom(); i+=5) {
      os << indent;