C++ PDB::remarks方法代码示例

int main(int argc, char *argv[]) {
  string hdr = invocationHeader(argc, argv);
  opts::BasicOptions* bopts = new opts::BasicOptions;
  opts::ModelWithCoords* mopts = new opts::ModelWithCoords;
  ToolOptions* topts = new ToolOptions;
  opts::AggregateOptions options;
  options.add(bopts).add(mopts).add(topts);

  if (!options.parse(argc, argv))
    exit(-1);

  AtomicGroup model = mopts->model;
  
  if (topts->reimage) {
    if (!model.isPeriodic()) {
      cerr << "WARNING- Reimaging requested, but the model has no periodic box information\n";
    } else {
      if (!topts->bonds_name.empty()) {
        AtomicGroup bonds = createSystem(topts->bonds_name);
        copyBonds(model, bonds);
      }

      if (!model.hasBonds()) {
        cerr << "WARNING- The model has no connectivity.  Assigning bonds based on distance.\n";
        model.findBonds();
      }
    }
  }

  AtomicGroup center_mol = selectAtoms(model, topts->center_sel);
  GCoord center = center_mol.centroid();

  AtomicGroup apply_mol = selectAtoms(model, topts->apply_sel);
  GCoord offset = topts->translate - center;
  if (topts->center_xy)
    offset.z() = 0.0;
  
  for (AtomicGroup::iterator atom = apply_mol.begin(); atom != apply_mol.end(); ++atom)
    (*atom)->coords() += offset;

  if (topts->reimage && model.isPeriodic()) {
    vGroup molecules = model.splitByMolecule();
    vGroup segments = model.splitByUniqueSegid();
      
    for (vGroup::iterator seg = segments.begin(); seg != segments.end(); ++seg)
      seg->reimage();

    for (vGroup::iterator mol =molecules.begin(); mol != molecules.end(); ++mol)
      mol->reimage();
  }

  AtomicGroup write_mol = selectAtoms(model, topts->write_sel);
  PDB pdb = PDB::fromAtomicGroup(write_mol);
  pdb.remarks().add(hdr);
  cout << pdb;
}

int main(int argc, char *argv[]) {

  string hdr = invocationHeader(argc, argv);
  opts::BasicOptions* bopts = new opts::BasicOptions(fullHelpMessage());
  ToolOptions* topts = new ToolOptions;
  opts::RequiredArguments* ropts = new opts::RequiredArguments;
  ropts->addArgument("matrix", "matrix-file");

  opts::AggregateOptions options;
  options.add(bopts).add(topts).add(ropts);
  if (!options.parse(argc, argv))
    exit(-1);

  string matrix_name = ropts->value("matrix");

  RealMatrix A;
  readAsciiMatrix(matrix_name, A);
  uint m = A.rows();

  if (rows.empty())
    for (uint i=0; i<m; ++i)
      rows.push_back(i);

  uint resid = 1;
  uint total_chunks = rows.size() / (chunksize ? chunksize : rows.size());
  uint chunk = 1;
  AtomicGroup model;

  // Warn about pymol...
  if (rows.size() >= 100000 && bonds)
    cerr << "Warning- PyMol may not correctly handle the CONECT records generated by this tool.\n";
  

  for (uint atomid = 0; atomid < rows.size(); ++atomid, ++resid) {
    if (chunksize && resid > chunksize) {
      if (bonds) {
        for (uint i=atomid - chunksize; i<atomid-1; ++i)
          model[i]->addBond(model[i+1]);
      }

      resid = 1;
      ++chunk;
    }

    GCoord c(scales[0] * A(rows[atomid], columns[0]),
             scales[1] * A(rows[atomid], columns[1]),
             scales[2] * A(rows[atomid], columns[2]));
    pAtom pa(new Atom(atomid+1, atom_name, c));
    pa->resid(resid);
    pa->resname(residue_name);

    ostringstream segstream;
    segstream << boost::format(segid_fmt) % chunk;
    pa->segid(segstream.str());

    double b;
    double q = 0.0;
    
    if (chunksize) {
      b = (100.0 * (resid-1)) / chunksize;
      q = static_cast<double>(chunk-1) / total_chunks;
    } else
      b = 100.0 * atomid / rows.size();

    pa->bfactor(b);
    pa->occupancy(q);

    model.append(pa);
  }
  
  if (bonds) {
      if (chunksize && resid > 1)
	  for (uint i=rows.size() - chunksize; i<rows.size()-1; ++i)
	      model[i]->addBond(model[i+1]);
      else if (!chunksize)
	  for (uint i=0; i<rows.size()-1; ++i)
	      model[i]->addBond(model[i+1]);
  }
  
  PDB pdb = PDB::fromAtomicGroup(model);
  pdb.remarks().add(hdr);
  cout << pdb;
}

int main(int argc, char *argv[]) {

  string hdr = invocationHeader(argc, argv);
  cout << "# " << hdr << endl;
  parseArgs(argc, argv);

  DoubleMatrix eigvals;
  readAsciiMatrix(eigvals_name, eigvals);

  DoubleMatrix eigvecs;
  readAsciiMatrix(eigvecs_name, eigvecs);

  if (modes.empty())
    for (uint i=0; i<eigvals.rows(); ++i)
      modes.push_back(i);

  uint n = modes.size();
  uint m = eigvecs.rows();

  // V = m x n matrix of eigenvectors
  DoubleMatrix V(m, n);
  for (uint i=0; i<n; ++i)
    for (uint j=0; j<m; ++j)
      V(j, i) = eigvecs(j, modes[i]);

  // Make a copy and scale by the appropriate eigenvalue,
  // remembering that eigenvalues are inverted for ENM,
  // or squaring and not inverting in the case of PCA
  DoubleMatrix VS = V.copy();
  for (uint i=0; i<n; ++i) {
    double e = eigvals[modes[i]];
    double s = (pca_input) ? (scale * e * e): (scale / e);
    for (uint j=0; j<m; ++j)
      VS(j, i) *= s;
  }
  
  // U = Covariance matrix
  DoubleMatrix U = loos::Math::MMMultiply(VS,V,false,true);

  // B-factors come from trace of diagonal 3x3 superblocks
  vector<double> B;
  double prefactor = 8.0 *  M_PI * M_PI / 3.0;
  for (uint i=0; i<m; i += 3) {
    double b = prefactor * (U(i,i) + U(i+1, i+1) + U(i+2, i+2));
    B.push_back(b);
    cout << boost::format("%-8d %g\n") % (i/3) % b;
  }

  if (!pdb_name.empty()) {
    AtomicGroup model = createSystem(pdb_name);
    AtomicGroup subset = selectAtoms(model, selection);

    if ((unsigned int)subset.size() != B.size()) {
      cerr << boost::format("Error- selection has %d atoms, but %d were expected.\n") % subset.size() % B.size();
      exit(-10);
    }

    for (uint i=0; i<B.size(); ++i)
      subset[i]->bfactor(B[i]);

    PDB pdb = PDB::fromAtomicGroup(model);
    pdb.remarks().add(hdr);
    ofstream ofs(out_name.c_str());
    ofs << pdb;
  }

}