LOOS  v2.3.2
vsa.cpp
1 /*
2  vsa
3 
4  Usage:
5  vsa [options] subset environment model output_prefix
6 
7  See:
8  Woodcock et al, J Chem Phys (2008) 129:214109
9  Haffner & Zheng, J Chem Phys (2009) 130:194111
10 
11 */
12 
13 
14 /*
15  This file is part of LOOS.
16 
17  LOOS (Lightweight Object-Oriented Structure library)
18  Copyright (c) 2009 Tod D. Romo
19  Department of Biochemistry and Biophysics
20  School of Medicine & Dentistry, University of Rochester
21 
22  This package (LOOS) is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation under version 3 of the License.
25 
26  This package is distributed in the hope that it will be useful,
27  but WITHOUT ANY WARRANTY; without even the implied warranty of
28  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29  GNU General Public License for more details.
30 
31  You should have received a copy of the GNU General Public License
32  along with this program. If not, see <http://www.gnu.org/licenses/>.
33 
34 */
35 
36 
37 
38 #include <loos.hpp>
39 
40 #include <limits>
41 
42 #include "hessian.hpp"
43 #include "enm-lib.hpp"
44 #include "vsa-lib.hpp"
45 
46 
47 using namespace std;
48 using namespace loos;
49 using namespace ENM;
50 
51 
52 namespace opts = loos::OptionsFramework;
53 namespace po = loos::OptionsFramework::po;
54 
55 
56 // Globals...
57 double normalization = 1.0;
58 double threshold = 1e-10;
59 
60 string hdr;
61 string subsystem_selection, environment_selection, model_name, prefix;
62 int verbosity = 0;
63 bool debug = false;
64 bool occupancies_are_masses;
65 string psf_file;
66 
67 string spring_desc;
68 bool nomass;
69 
70 
71 string fullHelpMessage() {
72 
73  string s =
74  "\n"
75  "SYNOPSIS\n"
76  "\n"
77  "Compute the vibrational subsystem analysis version of the\n"
78  "anisotropic network model. (See Woodcock et al.)\n"
79  "\n"
80  "DESCRIPTION\n"
81  "\n"
82  "Computes the VSA network model given a subsystem and an\n"
83  "environment selection. The output consists of several different\n"
84  "ASCII formatted matrices (that can be read by Matlab/Octave) and\n"
85  "depends on whether or not masses are included in the\n"
86  "calculation. If debugging is turned on (--debug), then the\n"
87  "intermediate matrices are written out:\n"
88  "\tfoo_H.asc - Composite Hessian\n"
89  "\tfoo_Hss.asc - Subsystem Hessian\n"
90  "\tfoo_Hee.asc - Environment Hessian\n"
91  "\tfoo_Hse.asc - Subsystem-Environment Hessian\n"
92  "\tfoo_Heei.asc - Inverted Environment Hessian\n"
93  "\tfoo_Hssp.asc - Effective Subsystem Hessian\n"
94  "\tfoo_Ms.asc - Subsystem mass (optional)\n"
95  "\tfoo_Me.asc - Environment mass (optional)\n"
96  "\tfoo_Msp.asc - Effective subsystem mass (optional)\n"
97  "\tfoo_R.asc - Cholesky decomposition of Msp (optional)\n"
98  "\n\n"
99  "* Unit Subsystem Mass, Zero Environment Mass *\n\n"
100  "Here, the effective subsystem Hessian is created and a Singular\n"
101  "Value Decomposition used to solve the eigenproblem:\n"
102  "\tfoo_U.asc = Subsystem eigenvectors\n"
103  "\tfoo_s.asc = Subsystem eigenvalues\n"
104  "\n\n"
105  "* Subsystem and Environment Mass *\n\n"
106  "The generalized eigenvalue problem is solved creating the\n"
107  "following matrices:\n"
108  "\tfoo_s.asc = Subsystem eigenvalues (mass-weighted)\n"
109  "\tfoo_U.asc = Subsystem eigenvectors (mass-weighted)\n"
110  "\n\n"
111  "* Spring Constant Control *\n\n"
112  "The spring constant used is controlled by the --spring option.\n"
113  "If only the name for the spring function is given, then the default\n"
114  "parameters are used. Alternatively, the name may include a\n"
115  "comma-separated list of parameters to be passed to the spring\n"
116  "function, i.e. --spring=distance,15.0\n\n"
117  "Available spring functions:\n";
118 
119  vector<string> names = springNames();
120  for (vector<string>::iterator i = names.begin(); i != names.end(); ++i)
121  s = s + "\t" + *i + "\n";
122 
123  s +=
124  "\n\n"
125  "* Mass Control *\n\n"
126  "VSA, by default, assumes that masses will be present. These can\n"
127  "come from one of two sources. If \"--psf foo.psf\" is given,\n"
128  "then the masses will be assigned using the \"foo.psf\" file. This\n"
129  "assumes that the atoms are in the same order between the PSF file\n"
130  "and the structure file given on the command line. Alternatively,\n"
131  "the occupancy field of the PDB can be used with the\n"
132  "\"--occupancies 1\" option. See the psf-masses tool for one way to\n"
133  "copy masses into a PDB's occupancies.\n"
134  "\n"
135  "To disable masses (i.e. use unit masses for the subsystem and\n"
136  "zero masses for the environment), use the \"--nomass 1\" option.\n"
137  "\n\n"
138  "EXAMPLES \n\n"
139  "\n"
140  "vsa --occupancies 1 foo.pdb 'segid == \"TRAN\" && name == \"CA\"'\\\n"
141  " 'segid != \"TRAN\" && name == \"CA\"' foo_vsa\n"
142  "\tCompute the VSA for a transmembrane region based on segid with the\n"
143  "\tmasses stored in the occupancy field of the PDB. Here the enviroment\n"
144  "\tcontains all other CA's in the system.\n"
145  "\n"
146  "vsa --psf foo.psf foo.pdb \"`cat selection` && name == 'CA'\" \\\n"
147  " \"not (`cat selection`) && name == 'CA'\" foo_vsa\n"
148  "\tCompute the VSA for a transmembrane region where the selection\n"
149  "\tis stored in a file and masses are taken from a PSF file.\n"
150  "\n"
151  "vsa --nomass 1 foo.pdb 'name == \"CA\"' 'name =~ \"^(C|O|N)$\"' foo_vsa\n"
152  "\tCompute the mass-less VSA with CAs as the subsystem and all other\n"
153  "\tbackbone atoms as the environment.\n"
154  "\n"
155  "vsa --nomass 1 --spring hca foo.pdb 'name == \"CA\"' 'name =~ \"^(C|O|N)$\"' foo_vsa\n"
156  "\tThe same example as above, but using the HCA spring constants.\n";
157 
158  return(s);
159 
160 }
161 
162 
163 
164 
165 // @cond TOOLS_INTERNAL
166 class ToolOptions : public opts::OptionsPackage {
167 public:
168 
169  void addGeneric(po::options_description& o) {
170  o.add_options()
171  ("psf", po::value<string>(&psf_file), "Take masses from the specified PSF file")
172  ("debug", po::value<bool>(&debug)->default_value(false), "Turn on debugging (output intermediate matrices)")
173  ("occupancies", po::value<bool>(&occupancies_are_masses)->default_value(false), "Atom masses are stored in the PDB occupancy field")
174  ("nomass", po::value<bool>(&nomass)->default_value(false), "Disable mass as part of the VSA solution")
175  ("spring,S", po::value<string>(&spring_desc)->default_value("distance"), "Spring method and arguments");
176  }
177 
178  string print() const {
179  ostringstream oss;
180  oss << boost::format("psf='%s', debug=%d, occupancies=%d, nomass=%d, spring='%s'")
181  % psf_file
182  % debug
183  % occupancies_are_masses
184  % nomass
185  % spring_desc;
186  return(oss.str());
187  }
188 
189 
190 };
191 // @endcond
192 
193 
194 
195 
196 int main(int argc, char *argv[]) {
197  hdr = invocationHeader(argc, argv);
198 
199  // Build up options
200  opts::BasicOptions* bopts = new opts::BasicOptions(fullHelpMessage());
202  ToolOptions* topts = new ToolOptions;
204  ropts->addArgument("subsystem", "subsystem-selection");
205  ropts->addArgument("environment", "environment-selection");
206  ropts->addArgument("prefix", "output-prefix");
207 
208  opts::AggregateOptions options;
209  options.add(bopts).add(mopts).add(topts).add(ropts);
210  if (!options.parse(argc, argv))
211  exit(-1);
212 
213  // Extract values
214  AtomicGroup model = mopts->model;
215  verbosity = bopts->verbosity;
216  subsystem_selection = ropts->value("subsystem");
217  environment_selection = ropts->value("environment");
218  prefix = ropts->value("prefix");
219 
220  // Ugly way of handling multiple methods for getting masses into the equation...
221  if (verbosity > 0)
222  cerr << "Assigning masses...\n";
223 
224  // Get masess from somewhere, or rely on the defaults provided by
225  // Atom (i.e. should be 1)
226  if (! psf_file.empty())
227  massFromPSF(model, psf_file);
228  else if (occupancies_are_masses)
229  massFromOccupancy(model);
230  else if (!nomass)
231  cerr << "WARNING- using default masses\n";
232 
233 
234  // Partition the model for building the composite Hessian
235  AtomicGroup subsystem = selectAtoms(model, subsystem_selection);
236  AtomicGroup environment = selectAtoms(model, environment_selection);
237  AtomicGroup composite = subsystem + environment;
238 
239  if (verbosity > 1) {
240  cerr << "Subsystem size is " << subsystem.size() << endl;
241  cerr << "Environment size is " << environment.size() << endl;
242  }
243 
244  // Determine which kind of scaling to apply to the Hessian...
245  SpringFunction* spring = 0;
246  spring = springFactory(spring_desc);
247 
248  SuperBlock* blocker = new SuperBlock(spring, composite);
249 
250  VSA vsa(blocker, subsystem.size());
251  vsa.prefix(prefix);
252  vsa.meta(hdr);
253  vsa.debugging(debug);
254  vsa.verbosity(verbosity);
255 
256  if (!nomass) {
257  DoubleMatrix M = getMasses(composite);
258  vsa.setMasses(M);
259  }
260 
261  vsa.solve();
262 
263  writeAsciiMatrix(prefix + "_U.asc", vsa.eigenvectors(), hdr, false);
264  writeAsciiMatrix(prefix + "_s.asc", vsa.eigenvalues(), hdr, false);
265 
266  // Be good...
267  delete spring;
268  delete blocker;
269 
270 }
Provides simple way to add command-line arguments (required options)
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
Namespace for encapsulating options processing.
This class creates superblocks in a hessian.
Definition: hessian.hpp:52
void addArgument(const std::string &name, const std::string &description)
Add a required argument given a name (tag) and a description (currently unused)
STL namespace.
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
DoubleMatrix getMasses(const AtomicGroup &grp)
Build the 3n x 3n diagonal mass matrix for a group.
Definition: enm-lib.cpp:65
Options common to all tools (including –fullhelp)
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
std::string value(const std::string &s) const
Retrieve the value for an argument.
Interface for ENM spring functions.
std::ostream & writeAsciiMatrix(std::ostream &os, const Math::Matrix< T, P, S > &M, const std::string &meta, const Math::Range &start, const Math::Range &end, const bool trans=false, F fmt=F())
Write a submatrix to a stream.
SpringFunction * springFactory(const std::string &spring_desc)
Factory function for generating new SpringFunction instances based on a user string.
void massFromOccupancy(AtomicGroup &grp)
The masses are stored in the occupancy field of a PDB...
Definition: enm-lib.cpp:58
Namespace to encapsulate Elastic Network Model routines.
Definition: anm-lib.hpp:32
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
std::vector< std::string > springNames()
List of possible names for springFactory()
Request a model with coordinates.
Vibrational subsystem analysis ENM.
Definition: vsa-lib.hpp:59
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
void massFromPSF(AtomicGroup &grp, const string &name)
Copy the masses from a PSF onto a group.
Definition: enm-lib.cpp:51
Combines a set of OptionsPackages.
AggregateOptions & add(OptionsPackage *pack)
Add a pointer to an OptionsPackage that will be used for options.
void prefix(const std::string &s)
Filename prefix when we have to write something out.
Definition: enm-lib.hpp:85