LOOS  v2.3.2
traj_transform.cpp
1 /*
2  traj_transform.cpp
3 
4  (c) 2011 Tod D. Romo, Grossfield Lab
5  Department of Biochemistry
6  University of Rochster School of Medicine and Dentistry
7 
8 
9  C++ template for writing a tool that transforms a subset of a trajectory
10 */
11 
12 
13 /*
14  This file is part of LOOS.
15 
16  LOOS (Lightweight Object-Oriented Structure library)
17  Copyright (c) 2011 Tod D. Romo
18  Department of Biochemistry and Biophysics
19  School of Medicine & Dentistry, University of Rochester
20 
21  This package (LOOS) is free software: you can redistribute it and/or modify
22  it under the terms of the GNU General Public License as published by
23  the Free Software Foundation under version 3 of the License.
24 
25  This package is distributed in the hope that it will be useful,
26  but WITHOUT ANY WARRANTY; without even the implied warranty of
27  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28  GNU General Public License for more details.
29 
30  You should have received a copy of the GNU General Public License
31  along with this program. If not, see <http://www.gnu.org/licenses/>.
32 
33 */
34 
35 
36 
37 #include <loos.hpp>
38 
39 using namespace std;
40 using namespace loos;
41 
42 namespace opts = loos::OptionsFramework;
43 namespace po = loos::OptionsFramework::po;
44 
45 
46 
47 // ----------------------------------------------------------------
48 // ***EDIT***
49 // The following code is for implementing tool-specific
50 // options if the tool requires them. If not, this section can be
51 // deleted.
52 
53 double option1;
54 int option2;
55 
56 
57 // The following conditional prevents this class from appearing in the
58 // Doxygen documentation for LOOS:
59 //
60 // @cond TOOLS_INTERNAL
61 class ToolOptions : public opts::OptionsPackage {
62 public:
63 
64  // Change these options to reflect what your tool needs
65  void addGeneric(po::options_description& o) {
66  o.add_options()
67  ("option1", po::value<double>(&option1)->default_value(0.0), "Tool Option #1")
68  ("option2", po::value<int>(&option2)->default_value(42), "Tool option #2");
69  }
70 
71  // The print() function returns a string that describes what all the
72  // options are set to (for logging purposes)
73  string print() const {
74  ostringstream oss;
75  oss << boost::format("option1=%f, option2=%d") % option1 % option2;
76  return(oss.str());
77  }
78 
79 };
80 // @endcond
81 // ----------------------------------------------------------------
82 
83 
84 
85 
86 
87 int main(int argc, char *argv[]) {
88 
89  // Store the invocation information for logging later
90  string header = invocationHeader(argc, argv);
91 
92  // Build up the command-line options for this tool by instantiating
93  // the appropriate OptionsPackage objects...
94 
95  // Basic options should be used by all tools. It provides help,
96  // verbosity, and the ability to read options from a config file
98 
99  // Require an output prefix
101 
102  // This tool can operate on a subset of atoms. The BasicSelection
103  // object provides the "--selection" option.
105 
106  // The BasicTrajectory object handles specifying a trajectory as
107  // well as a "--skip" option that lets the tool skip the first
108  // number of frames (i.e. equilibration). It creates a pTraj object
109  // that is already primed for reading...
111 
112  // ***EDIT***
113  // Tool-specific options can be included here...
114  ToolOptions* topts = new ToolOptions;
115 
116  // ***EDIT***
117  // All of the OptionsPackages are combined via the AggregateOptions
118  // object. First instantiate it, then add the desired
119  // OptionsPackage objects. The order is important. We recommend
120  // you progress from general (Basic and Selection) to more specific
121  // (model) and finally the tool options.
122  opts::AggregateOptions options;
123  options.add(bopts).add(oopts).add(sopts).add(tropts).add(topts);
124 
125  // Parse the command-line. If an error occurred, help will already
126  // be displayed and it will return a FALSE value.
127  if (!options.parse(argc, argv))
128  exit(-1);
129 
130  // Pull the model from the options object (it will include coordinates)
131  AtomicGroup model = tropts->model;
132 
133  // Pull out the trajectory...
134  pTraj traj = tropts->trajectory;
135 
136  // Select the desired atoms to operate over...
137  AtomicGroup subset = selectAtoms(model, sopts->selection);
138 
139  // Extract the output prefix
140  string prefix = oopts->prefix;
141 
142  // Set-up the DCD writer
143  DCDWriter outdcd(prefix + ".dcd");
144 
145  // Now iterate over all frames in the trajectory (excluding the skip
146  // region).
147  // Track whether we're on the first frame or not, for generating a
148  // PDB that corresponds to this trajectory...
149  bool first_frame = true;
150  while (traj->readFrame()) {
151 
152  // Update the coordinates ONLY for the subset of atoms we're
153  // interested in...
154  traj->updateGroupCoords(subset);
155 
156  // ***EDIT***
157  // Perform some transformation here
158 
159  // Write out the frame to the DCD
160  outdcd.writeFrame(subset);
161 
162  // If this is the first frame, then also write it out as a PDB
163  if (first_frame) {
164  first_frame = false;
165 
166  // Note: the atomids must be renumbered to be sequential. This
167  // can also break connectivity, so we prune the connectivity
168  // first, then renumber. We also work with a *copy* of the
169  // subset so that we don't alter the one used for reading in the
170  // trajectory.
171 
172  AtomicGroup frame_copy = subset.copy();
173  frame_copy.pruneBonds();
174  frame_copy.renumber();
175 
176  // Now convert to a PDB
177  PDB pdb = PDB::fromAtomicGroup(frame_copy);
178  pdb.remarks().add(header);
179  string output_pdbname = prefix + ".pdb";
180  ofstream ofs(output_pdbname.c_str());
181  if (!ofs) {
182  cerr << "Error: unable to open output file for PDB\n";
183  exit(-10);
184  }
185  ofs << pdb;
186  }
187 
188  }
189 
190 }
Remarks & remarks(void)
Accessor for the remarks object...
Definition: pdb.cpp:547
Gets a string as prefix for output files (–prefix)
pTraj trajectory
The trajectory, primed by the –skip value (if specified)
void add(const std::string s)
Add a remark.
Definition: pdb_remarks.cpp:52
Namespace for encapsulating options processing.
Provides a single LOOS selection (–selection)
STL namespace.
AtomicGroup model
Model that describes the trajectory.
void renumber(const int start=1, const int stride=1)
Renumber the atomid's of the contained atoms...
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
Options common to all tools (including –fullhelp)
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
PDB reading/writing class.
Definition: pdb.hpp:69
Basic trajectory with a –skip option.
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
AtomicGroup copy(void) const
Creates a deep copy of this group.
Definition: AtomicGroup.cpp:56
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
void pruneBonds()
Attempt to prune connectivity (only retain bonds to atoms within this AtomicGroup) ...
Namespace for most things not already encapsulated within a class.
Combines a set of OptionsPackages.
AggregateOptions & add(OptionsPackage *pack)
Add a pointer to an OptionsPackage that will be used for options.