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