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