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 "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 : {
103 : public:
104 : static void registerKeywords( Keywords& keys );
105 : explicit PdbRenumber(const CLToolOptions& co );
106 : int main(FILE* in, FILE*out,Communicator& pc) override;
107 4 : std::string description()const override {
108 4 : return "Modify atom numbers in a PDB, possibly using hybrid-36 coding";
109 : }
110 : };
111 :
112 10426 : PLUMED_REGISTER_CLTOOL(PdbRenumber,"pdbrenumber")
113 :
114 3473 : void PdbRenumber::registerKeywords( Keywords& keys ) {
115 3473 : CLTool::registerKeywords( keys );
116 6946 : keys.add("compulsory","--ipdb","specify the name of the input PDB file");
117 6946 : keys.add("compulsory","--opdb","specify the name of the output PDB file");
118 10419 : keys.add("optional","--firstatomnumber","specify the desired serial number of the first atom of the output file");
119 6946 : keys.add("optional","--atomnumbers","specify the desired serial numbers of the atoms of the output file using a separate list");
120 3473 : }
121 :
122 7 : PdbRenumber::PdbRenumber(const CLToolOptions& co ):
123 7 : CLTool(co)
124 : {
125 7 : inputdata=commandline;
126 7 : }
127 :
128 3 : int PdbRenumber::main(FILE* in, FILE*out,Communicator& pc) {
129 :
130 : std::string ipdb;
131 6 : parse("--ipdb",ipdb);
132 : std::string opdb;
133 3 : parse("--opdb",opdb);
134 :
135 3 : unsigned iat=0;
136 :
137 6 : parse("--firstatomnumber",iat);
138 :
139 : std::string atomnumbers;
140 6 : parse("--atomnumbers",atomnumbers);
141 :
142 3 : plumed_assert(ipdb.length()>0) << "please specify the input PDB with --ipdb";
143 3 : plumed_assert(opdb.length()>0) << "please specify the onput PDB with --opdb";
144 : std::fprintf(out," with input PDB: %s\n",ipdb.c_str());
145 : std::fprintf(out," with output PDB: %s\n",opdb.c_str());
146 :
147 : std::vector<unsigned> serials;
148 :
149 3 : if(atomnumbers.length()>0) {
150 1 : plumed_assert(iat==0) << "it is not possible to use both --atomnumbers and --firstatomnumber";
151 : std::fprintf(out," reading atom numbers from file %s\n",atomnumbers.c_str());
152 1 : IFile ifile;
153 1 : ifile.open(atomnumbers);
154 : std::string line;
155 7 : while(ifile.getline(line)) {
156 : int i;
157 6 : Tools::convert(line,i);
158 6 : serials.push_back(i);
159 : }
160 1 : } else {
161 2 : if(iat==0) iat=1;
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") iat=0;
190 0 : ofile << line << "\n";
191 : }
192 : }
193 :
194 3 : return 0;
195 3 : }
196 : }
197 :
198 : } // End of namespace
|
