hess.cc

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

    }
  Ref<SCBlockInfo> bi = new SCBlockInfo(total, nirrep, dims);
  for (i=0; i<nirrep; i++) {
      bi->set_subdim(i, symmbasis[i]->coldim());
    }
  RefSCDimension dsym = new SCDimension(bi);

  RefSCDimension bd3natom = new SCDimension(3*mol->natom());
  bd3natom->blocks()->set_subdim(0,d3natom);

  Ref<SCMatrixKit> symkit = new BlockedSCMatrixKit(kit);
  RefSCMatrix result(dsym, bd3natom, symkit);
  BlockedSCMatrix *bresult = dynamic_cast<BlockedSCMatrix*>(result.pointer());

  // put the symmetric basis in the result matrix
  for (i=0; i<nirrep; i++) {
      if (dims[i]>0) bresult->block(i).assign(symmbasis[i].t());
    }

  delete[] symmbasis;

  for (i=0; i<natom; i++) delete[] atom_map[i];
  delete[] atom_map;

  delete[] dims;
  return result;
}

void
MolecularHessian::set_energy(const Ref<MolecularEnergy> &)
{
}

MolecularEnergy*
MolecularHessian::energy() const
{
  return 0;
}

void
MolecularHessian::write_cartesian_hessian(const char *filename,
                                          const Ref<Molecule> &mol,
                                          const RefSymmSCMatrix &hess)
{
  int ntri = (3*mol->natom()*(3*mol->natom()+1))/2;
  double *hessv = new double[ntri];
  hess->convert(hessv);
  if (MessageGrp::get_default_messagegrp()->me() == 0) {
      int i,j;
      ofstream out(filename);
      // file format is version text 1
      out << "Hessian VT1" << endl;
      out << mol->natom() << " atoms" << endl;
      for (i=0; i<mol->natom(); i++) {
          out << scprintf("%2d % 15.12f % 15.12f % 15.12f",
                          mol->Z(i), mol->r(i,0), mol->r(i,1), mol->r(i,2))
              << endl;
          
        }
      const int nrow = 5;
      for (i=0; i<ntri; ) {
          for (j=0; j<nrow && i<ntri; j++,i++) {
              if (j>0) out << " ";
              out << scprintf("% 15.12f", hessv[i]);
            }
          out << endl;
        }
      out << "End Hessian" << endl;
    }
  delete[] hessv;
}

void
MolecularHessian::read_cartesian_hessian(const char *filename,
                                         const Ref<Molecule> &mol,
                                         const RefSymmSCMatrix &hess)
{
  int ntri = (3*mol->natom()*(3*mol->natom()+1))/2;
  double *hessv = new double[ntri];
  Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
  if (grp->me() == 0) {
      int i;
      ifstream in(filename);
      const int nline = 100;
      char linebuf[nline];
      in.getline(linebuf, nline);
      if (strcmp(linebuf,"Hessian VT1")) {
          ExEnv::errn() << "MolecularHessian: not given a hessian file" << endl;
          abort();
        }
      int natom;
      in >> natom;
      if (natom != mol->natom()) {
          ExEnv::errn() << "MolecularHessian: wrong number of atoms in hessianfile"
               << endl;
          abort();
        }
      in.getline(linebuf,nline);
      //ExEnv::outn() << "READ: should be atoms: " << linebuf << endl;
      for (i=0; i<mol->natom(); i++) {
          int Z;
          double x, y, z;
          in >> Z >> x >> y >> z;
          //ExEnv::outn() << "READ: " << Z << " " << x << " " << y << " " << z << endl;
        }
      for (i=0; i<ntri; i++) {
          in >> hessv[i];
          //ExEnv::outn() << "READ: hess[" << i << "] = " << hessv[i] << endl;
        }
      in.getline(linebuf, nline);
      //ExEnv::outn() << "READ: last line = " << linebuf << endl;
      if (strcmp(linebuf,"End Hessian")) {
          // try once more since there could be a left over new line
          in.getline(linebuf, nline);
          if (strcmp(linebuf,"End Hessian")) {
              //ExEnv::outn() << "READ: last line = " << linebuf << endl;
              ExEnv::errn() << "MolecularHessian: hessian file seems to be truncated"
                   << endl;
              abort();
            }
        }
    }
  grp->bcast(hessv,ntri);
  hess->assign(hessv);
  delete[] hessv;
}

/////////////////////////////////////////////////////////////////
// ReadMolecularHessian

static ClassDesc ReadMolecularHessian_cd(
  typeid(ReadMolecularHessian),"ReadMolecularHessian",1,"public MolecularHessian",
  0, create<ReadMolecularHessian>, create<ReadMolecularHessian>);

ReadMolecularHessian::ReadMolecularHessian(const Ref<KeyVal>&keyval):
  MolecularHessian(keyval)
{
  KeyValValueString default_filename(SCFormIO::fileext_to_filename(".hess"),
                                     KeyValValueString::Steal);
  filename_ = keyval->pcharvalue("filename", default_filename);
}

ReadMolecularHessian::ReadMolecularHessian(StateIn&s):
  SavableState(s),
  MolecularHessian(s)
{
  s.getstring(filename_);
}

ReadMolecularHessian::~ReadMolecularHessian()
{
  delete[] filename_;
}

void
ReadMolecularHessian::save_data_state(StateOut&s)
{
  MolecularHessian::save_data_state(s);
  s.putstring(filename_);
}

RefSymmSCMatrix
ReadMolecularHessian::cartesian_hessian()
{
  RefSymmSCMatrix hess = matrixkit()->symmmatrix(d3natom());
  read_cartesian_hessian(filename_, mol_, hess);
  return hess;
}

/////////////////////////////////////////////////////////////////
// GuessMolecularHessian

static ClassDesc GuessMolecularHessian_cd(
  typeid(GuessMolecularHessian),"GuessMolecularHessian",1,"public MolecularHessian",
  0, create<GuessMolecularHessian>, create<GuessMolecularHessian>);

GuessMolecularHessian::GuessMolecularHessian(const Ref<KeyVal>&keyval):
  MolecularHessian(keyval)
{
  coor_ << keyval->describedclassvalue("coor");
  if (mol_.null()) mol_ = coor_->molecule();
}

GuessMolecularHessian::GuessMolecularHessian(StateIn&s):
  SavableState(s),
  MolecularHessian(s)
{
  coor_ << SavableState::restore_state(s);
}

GuessMolecularHessian::~GuessMolecularHessian()
{
}

void
GuessMolecularHessian::save_data_state(StateOut&s)
{
  MolecularHessian::save_data_state(s);
  SavableState::save_state(coor_.pointer(),s);
}

RefSymmSCMatrix
GuessMolecularHessian::cartesian_hessian()
{
  RefSymmSCMatrix hessian(coor_->dim(), coor_->matrixkit());
  coor_->guess_hessian(hessian);
  RefSymmSCMatrix xhessian(coor_->dim_natom3(), coor_->matrixkit());
  coor_->to_cartesian(xhessian,hessian);
  return xhessian;
}

/////////////////////////////////////////////////////////////////
// DiagMolecularHessian

static ClassDesc DiagMolecularHessian_cd(
  typeid(DiagMolecularHessian),"DiagMolecularHessian",1,"public MolecularHessian",
  0, create<DiagMolecularHessian>, create<DiagMolecularHessian>);

DiagMolecularHessian::DiagMolecularHessian(const Ref<KeyVal>&keyval):
  MolecularHessian(keyval)
{
  diag_ = keyval->doublevalue("coor",KeyValValuedouble(1.0));
}

DiagMolecularHessian::DiagMolecularHessian(StateIn&s):
  SavableState(s),
  MolecularHessian(s)
{
  s.get(diag_);
}

DiagMolecularHessian::~DiagMolecularHessian()
{
}

void
DiagMolecularHessian::save_data_state(StateOut&s)
{
  MolecularHessian::save_data_state(s);
  s.put(diag_);
}

RefSymmSCMatrix
DiagMolecularHessian::cartesian_hessian()
{
  RefSymmSCMatrix xhessian(d3natom(), matrixkit());
  xhessian->assign(0.0);
  xhessian->shift_diagonal(diag_);
  return xhessian;
}

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

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