LOOS  v2.3.2
ufidpick.cpp
1 /*
2  ufidpick
3 
4  Picks uniform fiducial structures for a structural histogram, a la Lyman &
5  Zuckerman, J Phys Chem B (2007) 111:12876-12882
6 
7 
8  Usage- ufidpick model trajectory selection output-name probability [seed]
9 
10 */
11 
12 
13 
14 /*
15 
16  This file is part of LOOS.
17 
18  LOOS (Lightweight Object-Oriented Structure library)
19  Copyright (c) 2010, Tod D. Romo
20  Department of Biochemistry and Biophysics
21  School of Medicine & Dentistry, University of Rochester
22 
23  This package (LOOS) is free software: you can redistribute it and/or modify
24  it under the terms of the GNU General Public License as published by
25  the Free Software Foundation under version 3 of the License.
26 
27  This package is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with this program. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 
37 
38 #include <loos.hpp>
39 
40 #include "fid-lib.hpp"
41 
42 using namespace std;
43 using namespace loos;
44 
45 
46 string fullHelpMessage(void) {
47  string msg =
48  "\n"
49  "SYNOPSIS\n"
50  "\tPick fiducial structures for a structural histogram\n"
51  "\n"
52  "DESCRIPTION\n"
53  "\n"
54  "\tThis tool picks structures from a trajectory for use as fiducials in a\n"
55  "structural histogram. They are picked using bins with a uniform probability. For\n"
56  "more details, see Lyman & Zuckerman, J Phys Chem B (2007) 111:12876-12882.\n"
57  "\n"
58  "EXAMPLES\n"
59  "\n"
60  "\tufidpick model.pdb simulation.dcd all 'name == \"CA\"' zuckerman 0.1 >ufidpick.log\n"
61  "This example uses bins with a probability of 0.1 (i.e. 10 bins), using only\n"
62  "the alpha-carbons. The output files include a log of what structures were \n"
63  "picked, stored in ufidpick.log, as well as a trajectory containing just the\n"
64  "fiducial structures in zuckerman.dcd and the corresponding model file in zuckerman.pdb\n"
65  "\n"
66  "SEE ALSO\n"
67  "\tassign_frames, hierarchy, effsize.pl, neff\n";
68 
69  return(msg);
70 }
71 
72 
73 int main(int argc, char *argv[]) {
74  string hdr = invocationHeader(argc, argv);
75 
76  if (argc < 7 || argc > 8) {
77  cerr << "Usage - " << argv[0] << " model trajectory range|all selection output-name cutoff [seed]\n";
78  cerr << fullHelpMessage();
79  exit(-1);
80  }
81 
82  int opti = 1;
83  AtomicGroup model = createSystem(argv[opti++]);
84  model.clearBonds();
85 
86  pTraj traj = createTrajectory(argv[opti++], model);
87  string range(argv[opti++]);
88  string selection(argv[opti++]);
89  AtomicGroup subset = selectAtoms(model, selection);
90  string outname(argv[opti++]);
91  double cutoff = strtod(argv[opti++], 0);
92 
93  uint seed;
94  if (opti < argc) {
95  seed = strtoul(argv[opti++], 0, 10);
96  rng_singleton().seed(static_cast<unsigned int>(seed)); // This addresses a bug in earlier BOOSTs
97  } else
98  seed = randomSeedRNG();
99 
100 
101  cout << "# " << hdr << endl;
102  cout << "# seed = " << seed << endl;
103 
104 
105  vector<uint> source_frames;
106  if (range == "all")
107  for (uint i=0; i<traj->nframes(); ++i)
108  source_frames.push_back(i);
109  else
110  source_frames = parseRangeList<uint>(range);
111 
112  vector<uint> frames = trimFrames(source_frames, cutoff);
113  if (frames.size() != source_frames.size())
114  cout << "# WARNING- truncated last " << source_frames.size() - frames.size() << " frames\n";
115 
116  boost::tuple<vecGroup, vecUint> result = pickFiducials(subset, traj, frames, cutoff);
117  cout << "# n\tref\n";
118  vecGroup fiducials = boost::get<0>(result);
119  vecUint id = boost::get<1>(result);
120 
121  for (uint i=0; i<id.size(); ++i)
122  cout << i << "\t" << id[i] << "\n";
123 
124 
125  DCDWriter dcd(outname + ".dcd", fiducials, hdr);
126 
127  PDB pdb = PDB::fromAtomicGroup(fiducials[0]);
128  pdb.remarks().add(hdr);
129  string fname = outname + ".pdb";
130  ofstream ofs(fname.c_str());
131  ofs << pdb;
132 
133 }
Remarks & remarks(void)
Accessor for the remarks object...
Definition: pdb.cpp:547
void add(const std::string s)
Add a remark.
Definition: pdb_remarks.cpp:52
STL namespace.
AtomicGroup createSystem(const std::string &filename)
Factory function for reading in structure files.
Definition: sfactories.cpp:115
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
PDB reading/writing class.
Definition: pdb.hpp:69
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
A very lightweight class for writing simple DCDs.
Definition: dcdwriter.hpp:48
uint randomSeedRNG(void)
Randomly seeds the RNG.
base_generator_type & rng_singleton(void)
Suite-wide random number generator singleton.
void clearBonds(void)
Remove any bonding information present in contained atoms.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.