Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "core/PlumedMain.h"
26 : #include "tools/PDB.h"
27 :
28 : #include <cstdio>
29 :
30 : namespace PLMD {
31 : namespace generic {
32 :
33 : namespace {
34 : // the values are passed to CONSTANT as text: they are written with all their digits,
35 : // so that they are read back exactly (Tools::convert would only keep 6 digits)
36 24808 : std::string fullPrecision(double v) {
37 : char buffer[32];
38 : std::snprintf(buffer,sizeof(buffer),"%.17g",v);
39 24808 : return buffer;
40 : }
41 : }
42 :
43 : //+PLUMEDOC COLVAR PDB2CONSTANT
44 : /*
45 : Create a constant value from a PDB input file
46 :
47 : This shortcut converts the contents of a PDB file to one or more [CONSTANT](CONSTANT.md) actions.
48 : Converting PDB files to Values in this way is useful because it means that when we implement methods, like those
49 : in the [refdist](module_refdist.md) or [mapping](module_mapping.md) modules or the [RMSD](RMSD.md) action, that calculate the distance between two
50 : configurations those two configurations are both stored in PLUMED values. The same code can thus be used to
51 : calculate the difference between the instantaneous configuration
52 : and a constant reference configuration that was read from a file or between two reference configuration.
53 :
54 : The following example illustrates how this action can be used to read a set of reference atomic positions
55 :
56 : ```plumed
57 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
58 : ref: PDB2CONSTANT REFERENCE=regtest/basic/rt19/test0.pdb
59 : ```
60 :
61 : You can see how the reference positions are converted to [CONSTANT](CONSTANT.md) action that outputs a vector by expanding the shortcut.
62 :
63 : You can also use this command to read in multiple reference positions as illustrated below:
64 :
65 : ```plumed
66 : #SETTINGS INPUTFILES=regtest/mapping/rt-pathtools-3/all.pdb
67 : ref: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pathtools-3/all.pdb
68 : ```
69 :
70 : The [CONSTANT](CONSTANT.md) that is created by this action is a matrix. Each row of the output matrix contains one set of reference positions.
71 : Notice also that if you have a PDB input which contains multiple reference configurations you can create a vector constant by using the `NUMBER`
72 : keyword to specify the particular configuration that you would like to read in as shown below:
73 :
74 : ```plumed
75 : #SETTINGS INPUTFILES=regtest/mapping/rt-pathtools-3/all.pdb
76 : ref: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pathtools-3/all.pdb NUMBER=4
77 : ```
78 :
79 : The input above will reads in the fourth configuration in the input PDB file only.
80 :
81 : ## The PDB file format
82 :
83 : PLUMED uses the PDB file format here and in several other places
84 :
85 : - To read the molecular structure ([MOLINFO](MOLINFO.md)).
86 : - To read reference conformations ([RMSD](RMSD.md), but also in methods such as [FIT_TO_TEMPLATE](FIT_TO_TEMPLATE.md), etc).
87 :
88 : The implemented PDB reader expects a file formatted correctly according to the
89 : [PDB standard](http://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html).
90 : In particular, the following columns are read from ATOM records
91 :
92 : ````
93 : columns | content
94 : 1-6 | record name (ATOM or HETATM)
95 : 7-11 | serial number of the atom (starting from 1)
96 : 13-16 | atom name
97 : 18-20 | residue name
98 : 22 | chain id
99 : 23-26 | residue number
100 : 31-38 | x coordinate
101 : 39-46 | y coordinate
102 : 47-54 | z coordinate
103 : 55-60 | occupancy
104 : 61-66 | beta factor
105 : ````
106 :
107 : The PLUMED parser is slightly more permissive than the official PDB format
108 : in the fact that the format of real numbers is not fixed. In other words,
109 : any real number that can be parsed is OK and the dot can be placed anywhere. However,
110 : __columns are interpret strictly__. A sample PDB should look like the following
111 :
112 : ````
113 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
114 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
115 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
116 : ````
117 :
118 : Notice that serial numbers need not to be consecutive. In the three-line example above,
119 : only the coordinates of three atoms are provided. This is perfectly legal and indicates to PLUMED
120 : that information about these atoms only is available. This could be both for structural
121 : information in [MOLINFO](MOLINFO.md), where the other atoms would have no name assigned, and for
122 : reference structures used in [RMSD](RMSD.md), where only the provided atoms would be used to compute RMSD.
123 :
124 : ## Including arguments in PDB files
125 :
126 : If you wish to specify reference values for PLUMED Values in the REMARKS of a PLUMED input file like this:
127 :
128 : ````
129 : REMARK t1=-4.3345
130 : REMARK t2=3.4725
131 : END
132 : ````
133 :
134 : You can read in these reference values by using the PDB2CONSTANT command as follows:
135 :
136 : ```plumed
137 : #SETTINGS INPUTFILES=regtest/mapping/rt-pathtools-4/epath.pdb
138 : t1: TORSION ATOMS=1,2,3,4
139 : t2: TORSION ATOMS=5,6,7,8
140 : t1_ref: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pathtools-4/epath.pdb ARG=t1
141 : t2_ref: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pathtools-4/epath.pdb ARG=t2
142 : ```
143 :
144 : In this case the input must define values with the labels that are being read in from the reference file
145 : and separate PDB2CONSTANT commands are required for reading in `t1` and `t2`. Furthermore, because the
146 : input PDB file contains multiple frames vectors containing all the values for `t1` and `t2` are output from
147 : the constant commands that are created by the shortcuts in the above input. If you want to read only one of the
148 : configurations in the input PDB file you can use a pdb with a single frame or the `NUMBER` keyword described above.
149 :
150 : If, for any reason, you want to read data from a PDB file that is not a reference value for one of the values defined in
151 : your PLUMED input file you use the NOARGS flag as shown below:
152 :
153 : ```plumed
154 : #SETTINGS INPUTFILES=regtest/mapping/rt-pathtools-4/epath.pdb
155 : t1_ref: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pathtools-4/epath.pdb NOARGS ARG=t1
156 : t2_ref: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pathtools-4/epath.pdb NOARGS ARG=t2
157 : ```
158 :
159 : ## Occupancy and beta factors
160 :
161 : PLUMED also reads the occupancy and beta factors from the input PDB files. However, these columns of data are
162 : given a very special meaning.
163 : In cases where the PDB structure is used as a reference for an alignment (that's the case
164 : for instance in [RMSD](RMSD.md) and in [FIT_TO_TEMPLATE](FIT_TO_TEMPLATE.md)), the occupancy column is used
165 : to provide the weight of each atom in the alignment. In cases where, perhaps after alignment,
166 : the displacement between running coordinates and the provided PDB is computed, the beta factors
167 : are used as weight for the displacement.
168 : Since setting the weights to zero is the same as __not__ including an atom in the alignment or
169 : displacement calculation, the two following reference files would be equivalent when used in an [RMSD](RMSD.md)
170 : calculation. First file:
171 :
172 : ````
173 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
174 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
175 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 0.00 0.00
176 : ````
177 :
178 : Second file:
179 :
180 : ````
181 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
182 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
183 : ````
184 :
185 : However notice that many extra atoms with zero weight might slow down the calculation, so
186 : removing lines is better than setting their weights to zero.
187 : In addition, weights for alignment need not to be equivalent to weights for displacement.
188 : Starting with PLUMED 2.7, if all the weights are set to zero they will be normalized to be equal to the
189 : inverse of the number of involved atoms. This means that it will be possible to use files with
190 : the weight columns set to zero obtaining a meaningful result. In previous PLUMED versions,
191 : setting all weights to zero was resulting in an error instead.
192 :
193 : ## Systems with more than 100k atoms
194 :
195 : Notice that it very likely does not make any sense to compute the [RMSD](RMSD.md) or any other structural
196 : deviation __using__ many atoms. However, if the protein for which you want to compute [RMSD](RMSD.md)
197 : has atoms with large serial numbers (e.g. because it is located __after__ solvent in the sorted list of atoms)
198 : you might end up with troubles with the limitations of the PDB format. Indeed, since there are 5
199 : columns available for atom serial number, this number cannot be larger than 99999.
200 : In addition, providing [MOLINFO](MOLINFO.md) with names associated to atoms with a serial larger than 99999 would be impossible.
201 :
202 : Since PLUMED 2.4 we allow the [hybrid 36](http://cci.lbl.gov/hybrid_36/) format to be used to specify atom numbers.
203 : This format is not particularly widespread, but has the nice feature that it provides a one-to-one mapping
204 : between numbers up to approximately 80 millions and strings with 5 characters, plus it is backward compatible
205 : for numbers smaller than 100000. This is not true for notations like the hex notation exported by VMD.
206 : Using the hybrid 36 format, the ATOM records for atom ranging from 99997 to 100002 would read like these:
207 :
208 : ````
209 : ATOM 99997 Ar X 1 45.349 38.631 15.116 1.00 1.00
210 : ATOM 99998 Ar X 1 46.189 38.631 15.956 1.00 1.00
211 : ATOM 99999 Ar X 1 46.189 39.471 15.116 1.00 1.00
212 : ATOM A0000 Ar X 1 45.349 39.471 15.956 1.00 1.00
213 : ATOM A0000 Ar X 1 45.349 38.631 16.796 1.00 1.00
214 : ATOM A0001 Ar X 1 46.189 38.631 17.636 1.00 1.00
215 : ````
216 :
217 : There are tools that can be found to translate from integers to strings and back using hybrid 36 format
218 : (a simple python script can be found [here](https://sourceforge.net/p/cctbx/code/HEAD/tree/trunk/iotbx/pdb/hybrid_36.py)).
219 : In addition, as of PLUMED 2.5, we provide a command line tool that can be used to renumber atoms in a PDB file.
220 :
221 : */
222 : //+ENDPLUMEDOC
223 :
224 : class PDB2Constant : public ActionShortcut {
225 : public:
226 : static void registerKeywords(Keywords& keys);
227 : explicit PDB2Constant(const ActionOptions&);
228 : };
229 :
230 : PLUMED_REGISTER_ACTION(PDB2Constant,"PDB2CONSTANT")
231 :
232 184 : void PDB2Constant::registerKeywords(Keywords& keys) {
233 184 : ActionShortcut::registerKeywords( keys );
234 184 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure");
235 184 : keys.add("compulsory","NUMBER","0","if there are multiple structures in the pdb file you can specify that you want the RMSD from a specific structure by specifying its place in the file here. If NUMBER=0 then the RMSD from all structures are computed");
236 184 : keys.addFlag("NOARGS",false,"the arguments that are being read from the PDB file are not in the plumed input");
237 368 : keys.addInputKeyword("optional","ARG","scalar/vector","read this single argument from the input rather than the atomic structure");
238 368 : keys.setValueDescription("scalar/vector/matrix","a value that is constructed from the information in the PDB file");
239 184 : keys.needsAction("CONSTANT");
240 184 : }
241 :
242 122 : PDB2Constant::PDB2Constant(const ActionOptions& ao):
243 : Action(ao),
244 122 : ActionShortcut(ao) {
245 : std::string input;
246 122 : parse("REFERENCE",input);
247 : unsigned frame;
248 122 : parse("NUMBER",frame);
249 122 : bool noargs=false;
250 122 : std::vector<std::string> argn;
251 122 : parseVector("ARG",argn);
252 122 : std::vector<Value*> theargs;
253 122 : if( argn.size()>0 ) {
254 75 : parseFlag("NOARGS",noargs);
255 75 : if( !noargs ) {
256 73 : ActionWithArguments::interpretArgumentList( argn, plumed.getActionSet(), this, theargs );
257 2 : } else if( argn.size()>1 ) {
258 0 : error("can only read one argument at a time from input pdb file");
259 : } else {
260 2 : log.printf(" reading argument %s from file \n", argn[0].c_str() );
261 : }
262 : }
263 122 : if( theargs.size()>1 ) {
264 0 : error("can only read one argument at a time from input pdb file");
265 : }
266 :
267 122 : FILE* fp=std::fopen(input.c_str(),"r");
268 : bool do_read=true;
269 122 : std::vector<double> vals;
270 122 : if(!fp) {
271 0 : plumed_merror("could not open reference file " + input);
272 : }
273 122 : unsigned natoms=0, nframes=0;
274 :
275 2431 : while ( do_read ) {
276 2431 : PDB mypdb;
277 2431 : do_read=mypdb.readFromFilepointer(fp,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength());
278 2431 : if( !do_read && nframes>0 ) {
279 : break ;
280 : }
281 :
282 2309 : if( natoms==0 ) {
283 1385 : natoms = mypdb.getPositions().size();
284 924 : } else if( mypdb.getPositions().size()!=natoms ) {
285 0 : plumed_merror("mismatch between sizes of reference configurations");
286 : }
287 :
288 2309 : if( nframes+1==frame || frame==0 ) {
289 1831 : std::vector<double> align( mypdb.getOccupancy() );
290 : double asum=0;
291 12969 : for(unsigned i=0; i<align.size(); ++i) {
292 11138 : asum += align[i];
293 : }
294 1831 : if( asum>epsilon ) {
295 863 : double iasum = 1 / asum;
296 12001 : for(unsigned i=0; i<align.size(); ++i) {
297 11138 : align[i] *= iasum;
298 : }
299 968 : } else if( mypdb.size()>0 ) {
300 0 : double iasum = 1 / mypdb.size();
301 0 : for(unsigned i=0; i<align.size(); ++i) {
302 0 : align[i] = iasum;
303 : }
304 : }
305 1831 : Vector center;
306 : center.zero();
307 12969 : for(unsigned i=0; i<mypdb.getPositions().size(); ++i) {
308 11138 : center += align[i]*mypdb.getPositions()[i];
309 : }
310 :
311 1831 : if( theargs.size()==0 && argn.size()==0 ) {
312 2444 : for(unsigned j=0; j<3; ++j) {
313 25419 : for(unsigned i=0; i<mypdb.getPositions().size(); ++i) {
314 23586 : vals.push_back( mypdb.getPositions()[i][j] - center[j] );
315 : }
316 : }
317 1220 : } else if( noargs ) {
318 20 : std::vector<double> argvals( 1 );
319 20 : if( !mypdb.getArgumentValue(argn[0], argvals ) ) {
320 0 : error("argument " + argn[0] + " was not set in pdb input");
321 : }
322 20 : vals.push_back( argvals[0] );
323 20 : } else {
324 1200 : std::vector<double> argvals( theargs[0]->getNumberOfValues() );
325 1200 : if( !mypdb.getArgumentValue(theargs[0]->getName(), argvals ) ) {
326 0 : error("argument " + theargs[0]->getName() + " was not set in pdb input");
327 : }
328 2402 : for(unsigned i=0; i<argvals.size(); ++i) {
329 1202 : vals.push_back( argvals[i] );
330 : }
331 1200 : }
332 1831 : }
333 2309 : nframes++;
334 2431 : }
335 122 : if( frame>0 ) {
336 68 : nframes=1;
337 : }
338 122 : std::fclose(fp);
339 122 : plumed_assert( vals.size()>0 );
340 122 : std::string valstr = " VALUES=" + fullPrecision( vals[0] );
341 24808 : for(unsigned i=1; i<vals.size(); ++i) {
342 49372 : valstr += "," + fullPrecision( vals[i] );
343 : }
344 122 : if( vals.size()>nframes ) {
345 : std::string nc, nr;
346 48 : Tools::convert( nframes, nr );
347 48 : Tools::convert( vals.size()/nframes, nc );
348 144 : readInputLine( getShortcutLabel() + ": CONSTANT NROWS=" + nr + " NCOLS=" + nc + valstr );
349 : } else {
350 148 : readInputLine( getShortcutLabel() + ": CONSTANT" + valstr );
351 : }
352 366 : }
353 :
354 : }
355 : }
|
