Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2018-2020 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 "CLTool.h"
23 : #include "CLToolRegister.h"
24 : #include "tools/Tools.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/IFile.h"
27 : #include "tools/OFile.h"
28 : #include "tools/h36.h"
29 : #include <cstdio>
30 : #include <string>
31 : #include <vector>
32 : #include <array>
33 : #include <limits>
34 :
35 : using namespace std;
36 :
37 : namespace PLMD {
38 : namespace cltools {
39 :
40 : //+PLUMEDOC TOOLS pdbrenumber
41 : /*
42 : Modify atom numbers in a PDB, possibly using hybrid-36 coding.
43 :
44 : When reading a PDB files, PLUMED honors the serial number of each atom.
45 : This command can be used to process a PDB file changing the atom serial numbers.
46 : Notice that the resulting list might have gaps. It is however fundamental
47 : that atom numbers correspond to those used within the MD code.
48 : Importantly, if the serial number of an atom is greater than 99999, it is
49 : written in hybrid-36 notation (see \ref pdbreader).
50 : The main use of \ref pdbrenumber is thus that of producing files where atoms
51 : are numbered using hybrid-36 convention.
52 :
53 : The output PDB file is identical to the input PDB file, except for the atom number
54 : field.
55 : The rest of the line is written unchanged
56 : to the output file, even if it is incorrectly formatted. Residue numbers are not touched,
57 : and atom numbers in the input file are ignored.
58 :
59 :
60 : \par Examples
61 :
62 : By default, \ref pdbreader just sets the numbers as progressive starting from 1.
63 : For instance the following command:
64 : \verbatim
65 : > plumed pdbrenumber --ipdb input.pdb --opdb output.pdb
66 : \endverbatim
67 : will copy file `input.pdb` to `output.pdb` replacing all the serial atoms with
68 : increasing numbers starting from one. Atoms that have an index that is greater than 99999 will be written
69 : in the output PDB file in hybrid-36 code.
70 :
71 : It is possible to set a different serial number for the first atom, letting the
72 : following ones grow by one at each line. Here for instance the first atom
73 : will be assigned serial 1000, the second serial 1001, etc:
74 : \verbatim
75 : > plumed pdbrenumber --ipdb input.pdb --opdb output.pdb --firstatomnumber 1000
76 : \endverbatim
77 : If the first atom number is >99999, it should be given as a decimal number (not in hybrid-36 code).
78 : However, numbers >99999 in the output PDB file will be written in hybrid-36 code.
79 :
80 : As an alternative, one can provide a list of atoms as one per line in an auxiliary file.
81 : \verbatim
82 : > plumed pdbrenumber --ipdb input.pdb --opdb output.pdb --atomnumbers list.txt
83 : \endverbatim
84 : The `list.txt` file might be something like this
85 : \verbatim
86 : 120000
87 : 120001
88 : 120002
89 : 1
90 : 2
91 : 3
92 : \endverbatim
93 : Numbers >99999 in the list should be provided as decimal numbers (not in hybrid-36 code).
94 : However, numbers >99999 in the output PDB file will be written in hybrid-36 code.
95 : Notice that there should be at least enough lines in `list.txt` as many atoms in the PDB file.
96 : Additional lines in `list.txt` will just be ignored.
97 :
98 :
99 : */
100 : //+ENDPLUMEDOC
101 :
102 3 : class PdbRenumber:
103 : public CLTool
104 : {
105 : public:
106 : static void registerKeywords( Keywords& keys );
107 : explicit PdbRenumber(const CLToolOptions& co );
108 : int main(FILE* in, FILE*out,Communicator& pc);
109 0 : string description()const {
110 0 : return "Modify atom numbers in a PDB, possibly using hybrid-36 coding";
111 : }
112 : };
113 :
114 7362 : PLUMED_REGISTER_CLTOOL(PdbRenumber,"pdbrenumber")
115 :
116 1839 : void PdbRenumber::registerKeywords( Keywords& keys ) {
117 1839 : CLTool::registerKeywords( keys );
118 7356 : keys.add("compulsory","--ipdb","specify the name of the input PDB file");
119 7356 : keys.add("compulsory","--opdb","specify the name of the output PDB file");
120 7356 : keys.add("optional","--firstatomnumber","specify the desired serial number of the first atom of the output file");
121 7356 : keys.add("optional","--atomnumbers","specify the desired serial numbers of the atoms of the output file using a separate list");
122 1839 : }
123 :
124 3 : PdbRenumber::PdbRenumber(const CLToolOptions& co ):
125 3 : CLTool(co)
126 : {
127 3 : inputdata=commandline;
128 3 : }
129 :
130 3 : int PdbRenumber::main(FILE* in, FILE*out,Communicator& pc) {
131 :
132 : std::string ipdb;
133 6 : parse("--ipdb",ipdb);
134 : std::string opdb;
135 6 : parse("--opdb",opdb);
136 :
137 3 : unsigned iat=0;
138 :
139 6 : parse("--firstatomnumber",iat);
140 :
141 : std::string atomnumbers;
142 6 : parse("--atomnumbers",atomnumbers);
143 :
144 3 : plumed_assert(ipdb.length()>0) << "please specify the input PDB with --ipdb";
145 3 : plumed_assert(opdb.length()>0) << "please specify the onput PDB with --opdb";
146 : fprintf(out," with input PDB: %s\n",ipdb.c_str());
147 : fprintf(out," with output PDB: %s\n",opdb.c_str());
148 :
149 : std::vector<unsigned> serials;
150 :
151 3 : if(atomnumbers.length()>0) {
152 1 : plumed_assert(iat==0) << "it is not possible to use both --atomnumbers and --firstatomnumber";
153 : fprintf(out," reading atom numbers from file %s\n",atomnumbers.c_str());
154 2 : IFile ifile;
155 1 : ifile.open(atomnumbers);
156 : std::string line;
157 13 : while(ifile.getline(line)) {
158 : int i;
159 6 : Tools::convert(line,i);
160 12 : serials.push_back(i);
161 : }
162 : } else {
163 2 : if(iat==0) iat=1;
164 2 : fprintf(out," with atoms starting from %u\n",iat);
165 : }
166 :
167 6 : IFile ifile;
168 3 : ifile.open(ipdb);
169 :
170 6 : OFile ofile;
171 3 : ofile.open(opdb);
172 :
173 : std::string line;
174 200221 : while(ifile.getline(line)) {
175 100109 : std::string record=line.substr(0,6);
176 100109 : Tools::trim(record);
177 :
178 100109 : if(record=="ATOM" || record=="HETATM") {
179 : std::array<char,6> at;
180 100109 : unsigned ii=iat;
181 100109 : if(serials.size()>0) {
182 6 : plumed_assert(iat<serials.size()) << "there are more atoms in the PDB than serials in the file";
183 6 : ii=serials[iat];
184 : }
185 100109 : const char* msg = h36::hy36encode(5,ii,&at[0]);
186 100109 : plumed_assert(msg==nullptr) << msg;
187 100109 : at[5]=0;
188 300327 : ofile << line.substr(0,6) << &at[0] << line.substr(11) << "\n";
189 100109 : iat++;
190 : } else {
191 0 : if(record=="END" || record=="ENDMDL") iat=0;
192 0 : ofile << line << "\n";
193 : }
194 : }
195 :
196 3 : return 0;
197 : }
198 : }
199 :
200 5517 : } // End of namespace
|
