LOOS  v2.3.2
avgconv.cpp
1 /*
2  avgconv
3 
4  Usage- avgconv model traj selection range
5 */
6 
7 
8 
9 /*
10 
11  This file is part of LOOS.
12 
13  LOOS (Lightweight Object-Oriented Structure library)
14  Copyright (c) 2011, Tod D. Romo
15  Department of Biochemistry and Biophysics
16  School of Medicine & Dentistry, University of Rochester
17 
18  This package (LOOS) is free software: you can redistribute it and/or modify
19  it under the terms of the GNU General Public License as published by
20  the Free Software Foundation under version 3 of the License.
21 
22  This package is distributed in the hope that it will be useful,
23  but WITHOUT ANY WARRANTY; without even the implied warranty of
24  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  GNU General Public License for more details.
26 
27  You should have received a copy of the GNU General Public License
28  along with this program. If not, see <http://www.gnu.org/licenses/>.
29 */
30 
31 
32 
33 
34 
35 #include <loos.hpp>
36 #include <boost/format.hpp>
37 
38 using namespace std;
39 using namespace loos;
40 
41 
42 bool locally_optimal = false;
43 
44 
45 
46 string fullHelpMessage(void) {
47  string msg =
48  "\n"
49  "SYNOPSIS\n"
50  "\tConvergence of the average structure\n"
51  "\n"
52  "DESCRIPTION\n"
53  "\n"
54  "\tThe convergence of the average structure from a trajectory is determined by first\n"
55  "dividing the trajectory into blocks. An average structure is computed for the i'th and\n"
56  "the i+1'th block. These two average structure are superimposed using a Kabsch alignment\n"
57  "algorithm. The RMSD is calculated. This is then repeated for all blocks.\n"
58  "\tInitially, the whole trajectory is optimally aligned first using an iterative alignment\n"
59  "process (described in Grossfield, et al. Proteins 67, 31–40 (2007)). Optionally,\n"
60  "each block may be optimally aligned independently by using the 'locally optimal' flag.\n"
61  "\n"
62  "EXAMPLES\n"
63  "\n"
64  "\tavgconv sim.psf traj.dcd 'segid == \"PE1\"' >avgconv.asc\n"
65  "This example uses automatic block-sizes for the subsamples and calculates the RMSD and\n"
66  "superpositions using all atoms from the PE1 segment. The output is placed in avgconv.asc\n"
67  "\n"
68  "\tavgconv sim.psf traj.dcd 'name == \"CA\"' 10:10:500 >avgconv.asc\n"
69  "This example uses all alpha-carbons and explicitly sets the block sizes to 10,\n"
70  "20, ..., 100\n"
71  "\n"
72  "\tavgconv sim.psf traj.dcd '!hydrogen' 10:10:500 1 >avgconv.asc\n"
73  "This example uses all non-hydrogen atoms with block sizes of 10, 20, 30, ..., 5000,\n"
74  "and the blocks are all iteratively aligned prior to computing the average.\n"
75  "SEE ALSO\n"
76  "\tblock_average, block_avgconv\n";
77 
78  return(msg);
79 }
80 
81 
82 AtomicGroup calcAverage(const vector<AtomicGroup>& ensemble, const uint size) {
83 
84  vector<AtomicGroup> subsample(size);
85  for (uint i=0; i<size; ++i)
86  subsample[i] = ensemble[i];
87 
88  if (locally_optimal)
89  (void)iterativeAlignment(subsample);
90 
91  AtomicGroup avg = averageStructure(subsample);
92  return(avg);
93 }
94 
95 
96 
97 int main(int argc, char *argv[]) {
98  if (argc < 4 || argc > 6) {
99  cerr << "Usage- avgconv model traj selection [range [1 = local optimal avg]]\n";
100  cerr << fullHelpMessage();
101  exit(-1);
102  }
103 
104  string hdr = invocationHeader(argc, argv);
105  cout << "# " << hdr << endl;
106  cout << "# n\trmsd\n";
107 
108  int k = 1;
109  AtomicGroup model = createSystem(argv[k++]);
110  pTraj traj = createTrajectory(argv[k++], model);
111  string sel = string(argv[k++]);
112 
113  vector<uint> blocks;
114  if (argc == 4) {
115 
116  uint step = traj->nframes() / 100;
117  if (step == 0) {
118  cerr << "Error- too few frames for auto-blocksizes. Please specify block sizes explicitly\n";
119  exit(-1);
120  }
121  for (uint i=step; i<traj->nframes(); i += step)
122  blocks.push_back(i);
123 
124  } else {
125 
126  blocks = parseRangeList<uint>(argv[k++]);
127  if (argc == k+1)
128  locally_optimal = (argv[k][0] == '1');
129  }
130 
131  AtomicGroup subset = selectAtoms(model, sel);
132  cout << boost::format("# Subset has %d atoms\n") % subset.size();
133  vector<AtomicGroup> ensemble;
134  readTrajectory(ensemble, subset, traj);
135  cout << boost::format("# Trajectory has %d frames\n") % ensemble.size();
136 
137  cout << boost::format("# Blocks = %d\n") % blocks.size();
138 
139  if (!locally_optimal) {
140  boost::tuple<vector<XForm>, greal, int> result = iterativeAlignment(ensemble);
141  cout << boost::format("# Iterative alignment converged to RMSD of %g with %d iterations\n") % boost::get<1>(result) % boost::get<2>(result);
142  }
143 
144  AtomicGroup preceding = calcAverage(ensemble, blocks[0]);
145  for (vector<uint>::const_iterator ci = blocks.begin()+1; ci != blocks.end(); ++ci) {
146  AtomicGroup avg = calcAverage(ensemble, *ci);
147  avg.alignOnto(preceding);
148  double rmsd = preceding.rmsd(avg);
149 
150  cout << boost::format("%d\t%f\n") % *ci % rmsd;
151 
152  preceding = avg;
153  }
154 
155 }
156 
def averageStructure(traj)
Python version of averageStructure (using loos.pyloos.Trajectory-like objects) The subset defined in ...
Definition: ensembles.py:14
def iterativeAlignment(ensemble, threshold=1e-8, maxiter=1000)
Iteratively align a loos.pyloos.Trajectory object (or a list of AtomicGroups)
Definition: alignment.py:12
STL namespace.
greal rmsd(const AtomicGroup &)
Compute the RMSD between two groups.
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
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
GMatrix alignOnto(const AtomicGroup &)
Superimposes the current group onto the passed group.
Definition: AG_linalg.cpp:196
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.