Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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/ActionWithArguments.h"
23 : #include "core/ActionWithValue.h"
24 : #include "core/ActionAtomistic.h"
25 : #include "core/ActionPilot.h"
26 : #include "core/ActionRegister.h"
27 : #include "core/PlumedMain.h"
28 : #include "tools/OFile.h"
29 : #include "tools/h36.h"
30 :
31 : namespace PLMD {
32 : namespace generic {
33 :
34 : //+PLUMEDOC PRINTANALYSIS DUMPPDB
35 : /*
36 : Output PDB file.
37 :
38 : This command can be used to output a PDB file that contains the result from a dimensionality reduction
39 : calculation. To understand how to use this command you are best off looking at the examples of dimensionality
40 : reduction calculations that are explained in the documentation of the actions in the [dimred](module_dimred.md) module.
41 :
42 : Please note that this command __cannot__ be used in place of [DUMPATOMS](DUMPATOMS.md) to output a pdb file rather than
43 : a gro or xyz file. This is an example where this command is used to output atomic positions:
44 :
45 : ```plumed
46 : ff: COLLECT_FRAMES ATOMS=1-22 STRIDE=1
47 : DUMPPDB ATOM_INDICES=1-22 ATOMS=ff_data FILE=pos.pdb
48 : ```
49 :
50 : From this example alone it is clear that this command works very differently to the [DUMPATOMS](DUMPATOMS.md) command.
51 :
52 : As indcated below, you can specify what should appear in the title line of your PDB file as well as the residue indices, occupancies
53 : and beta values for each of the collected atoms.
54 :
55 : ```plumed
56 : ff: COLLECT_FRAMES ATOMS=1-4 STRIDE=1 CLEAR=100
57 : DUMPPDB ...
58 : ATOM_INDICES=1-4 ATOMS=ff_data
59 : DESCRIPTION=title
60 : RESIDUE_INDICES=1,1,2,2
61 : OCCUPANCY=1,0.5,0.5,0.75
62 : BETA=0.5,1,0.2,0.6
63 : FILE=pos.pdb
64 : STRIDE=100
65 : ...
66 : ```
67 :
68 : Notice, furthermore, that we store 100 and output 100 frame blocks of data here from the trajectory. The input above will thus output
69 : multiple pdb files.
70 :
71 : Notice, furthermore, that you can output the positions of atoms along with the argument values that correspond to each
72 : set of atomic positions as follows:
73 :
74 : ```plumed
75 : # Calculate the distance between atoms 1 and 2
76 : d: DISTANCE ATOMS=1,2
77 : # Collect the distance between atoms 1 and 2
78 : cd: COLLECT ARG=d
79 : # Calculate the angle involving atoms 1, 2 and 3
80 : a: ANGLE ATOMS=1,2,3
81 : # Collect the angle involving atoms 1, 2 and 3
82 : ca: COLLECT ARG=a
83 : # Now collect the positions
84 : ff: COLLECT_FRAMES ATOMS=1-22 STRIDE=1
85 : # Output a PDB file that contains the collected atomic positions
86 : # and the collected distances and angles
87 : DUMPPDB ATOM_INDICES=1-22 ATOMS=ff_data ARG=cd,ca FILE=pos.pdb
88 : ```
89 :
90 : The outputted pdb file has the format described in the documentation for the [PDB2CONSTANT](PDB2CONSTANT.md) action.
91 : The names that are used for each of the arguments in the pdb file are the same as the labels of the values in the input above (`d` and `a` in the above case).
92 : If for any reason you want to use different names for each of the arguments in the PDB file you can use the `ARG_NAMES` keyword as
93 : shown below:
94 :
95 : ```plumed
96 : # Calculate the distance between atoms 1 and 2
97 : d: DISTANCE ATOMS=1,2
98 : # Collect the distance between atoms 1 and 2
99 : cd: COLLECT ARG=d
100 : # Calculate the angle involving atoms 1, 2 and 3
101 : a: ANGLE ATOMS=1,2,3
102 : # Collect the angle involving atoms 1, 2 and 3
103 : ca: COLLECT ARG=a
104 : # Now collect the positions
105 : ff: COLLECT_FRAMES ATOMS=1-22 STRIDE=1
106 : # Output a PDB file that contains the collected atomic positions
107 : # and the collected distances and angles
108 : DUMPPDB ATOM_INDICES=1-22 ATOMS=ff_data ARG=cd,ca ARG_NAMES=m,n FILE=pos.pdb
109 : ```
110 :
111 : The distance will be shown with the name `m` for the input above while the angle will have the name `n`.
112 :
113 : */
114 : //+ENDPLUMEDOC
115 :
116 : class DumpPDB :
117 : public ActionWithArguments,
118 : public ActionPilot {
119 : std::string fmt;
120 : std::string file;
121 : std::string description;
122 : std::vector<unsigned> pdb_resid_indices;
123 : std::vector<double> beta, occupancy;
124 : std::vector<std::string> argnames;
125 : std::vector<AtomNumber> pdb_atom_indices;
126 : void buildArgnames();
127 : void printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const ;
128 : public:
129 : static void registerKeywords( Keywords& keys );
130 : explicit DumpPDB(const ActionOptions&);
131 13 : void calculate() override {}
132 13 : void apply() override {}
133 : void update() override ;
134 : };
135 :
136 : PLUMED_REGISTER_ACTION(DumpPDB,"DUMPPDB")
137 :
138 25 : void DumpPDB::registerKeywords( Keywords& keys ) {
139 25 : Action::registerKeywords( keys );
140 25 : ActionPilot::registerKeywords( keys );
141 25 : ActionWithArguments::registerKeywords( keys );
142 50 : keys.addInputKeyword("optional","ARG","vector/matrix","the values that are being output in the PDB file");
143 25 : keys.add("optional","ATOMS","value containing positions of atoms that should be output");
144 25 : keys.add("compulsory","STRIDE","0","the frequency with which the atoms should be output");
145 25 : keys.add("hidden","FMT","the format that should be used to output real numbers in the title");
146 25 : keys.add("compulsory","FILE","the name of the file on which to output these quantities");
147 25 : keys.add("compulsory","OCCUPANCY","1.0","vector of values to output in the occupancy column of the pdb file");
148 25 : keys.add("compulsory","BETA","1.0","vector of values to output in the beta column of the pdb file");
149 25 : keys.add("optional","DESCRIPTION","the title to use for your PDB output");
150 25 : keys.add("optional","ATOM_INDICES","the indices of the atoms in your PDB output");
151 25 : keys.add("optional","RESIDUE_INDICES","the indices of the residues in your PDB output");
152 25 : keys.add("optional","ARG_NAMES","the names of the arguments that are being output");
153 25 : }
154 :
155 16 : DumpPDB::DumpPDB(const ActionOptions&ao):
156 : Action(ao),
157 : ActionWithArguments(ao),
158 16 : ActionPilot(ao) {
159 32 : parse("FILE",file);
160 16 : if(file.length()==0) {
161 0 : error("name out output file was not specified");
162 : }
163 16 : std::vector<std::string> atoms;
164 32 : parseVector("ATOMS",atoms);
165 16 : if( atoms.size()>0 ) {
166 9 : std::vector<Value*> atom_args;
167 9 : interpretArgumentList( atoms, plumed.getActionSet(), this, atom_args );
168 9 : if( atom_args.size()!=1 ) {
169 0 : error("only one action should appear in input to atom args");
170 : }
171 9 : std::vector<Value*> args( getArguments() );
172 9 : args.push_back( atom_args[0] );
173 9 : requestArguments( args );
174 9 : std::vector<std::string> indices;
175 9 : parseVector("ATOM_INDICES",indices);
176 9 : std::vector<Value*> xvals,yvals,zvals,masv,chargev;
177 9 : ActionAtomistic::getAtomValuesFromPlumedObject(plumed,xvals,yvals,zvals,masv,chargev);
178 9 : ActionAtomistic::interpretAtomList( indices, xvals, this, pdb_atom_indices );
179 9 : log.printf(" printing atoms : ");
180 167 : for(unsigned i=0; i<indices.size(); ++i) {
181 158 : log.printf("%d ", pdb_atom_indices[i].serial() );
182 : }
183 9 : log.printf("\n");
184 9 : if( pdb_atom_indices.size()!=atom_args[0]->getShape()[1]/3 ) {
185 0 : error("mismatch between size of matrix containing positions and vector of atom indices");
186 : }
187 9 : beta.resize( atom_args[0]->getShape()[1]/3 );
188 9 : occupancy.resize( atom_args[0]->getShape()[1]/3 );
189 9 : parseVector("OCCUPANCY", occupancy );
190 9 : parseVector("BETA", beta );
191 18 : parseVector("RESIDUE_INDICES",pdb_resid_indices);
192 9 : if( pdb_resid_indices.size()==0 ) {
193 5 : pdb_resid_indices.resize( pdb_atom_indices.size() );
194 111 : for(unsigned i=0; i<pdb_resid_indices.size(); ++i) {
195 106 : pdb_resid_indices[i] = pdb_atom_indices[i].serial();
196 : }
197 4 : } else if( pdb_resid_indices.size()!=pdb_atom_indices.size() ) {
198 0 : error("mismatch between number of atom indices provided and number of residue indices provided");
199 : }
200 9 : }
201 16 : log.printf(" printing configurations to PDB file to file named %s \n", file.c_str() );
202 16 : if( getNumberOfArguments()==0 ) {
203 0 : error("nothing to print: use ARG and/or ATOMS");
204 : }
205 32 : parseVector("ARG_NAMES",argnames);
206 16 : if( argnames.size()==0 ) {
207 15 : buildArgnames();
208 : }
209 16 : fmt="%f";
210 16 : parse("FMT",fmt);
211 16 : fmt=" "+fmt;
212 16 : if( getStride()==0 ) {
213 8 : setStride(0);
214 8 : log.printf(" printing pdb at end of calculation \n");
215 : }
216 :
217 32 : parse("DESCRIPTION",description);
218 16 : if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->hasDerivatives() ) {
219 0 : error("argument for printing of PDB should be vector or matrix");
220 : }
221 16 : unsigned nv=getPntrToArgument(0)->getShape()[0];
222 46 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
223 30 : if( getPntrToArgument(i)->getRank()==0 || getPntrToArgument(i)->hasDerivatives() ) {
224 0 : error("argument for printing of PDB should be vector or matrix");
225 : }
226 30 : if( getPntrToArgument(i)->getShape()[0]!=nv ) {
227 0 : error("for printing to pdb files all arguments must have same number of values");
228 : }
229 : }
230 16 : }
231 :
232 2845 : void DumpPDB::printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const {
233 2845 : if( rnum>999 ) {
234 0 : plumed_merror("atom number too large to be used as residue number");
235 : }
236 : std::array<char,6> at;
237 2845 : const char* msg = h36::hy36encode(5,anum,&at[0]);
238 2845 : plumed_assert(msg==nullptr) << msg;
239 2845 : at[5]=0;
240 2845 : double lunits = plumed.getUnits().getLength()/0.1;
241 2845 : opdbf.printf("ATOM %s X RES %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
242 : &at[0], rnum, lunits*pos[0], lunits*pos[1], lunits*pos[2], m, q );
243 2845 : }
244 :
245 15 : void DumpPDB::buildArgnames() {
246 15 : unsigned nargs = getNumberOfArguments();
247 15 : if( pdb_atom_indices.size()>0 ) {
248 9 : nargs = nargs - 1;
249 : }
250 :
251 15 : argnames.resize(0);
252 15 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
253 : if( getPntrToArgument(0)->getRank()==2 ) {
254 : nvals = getPntrToArgument(0)->getShape()[0];
255 : }
256 35 : for(unsigned i=0; i<nargs; ++i) {
257 20 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
258 0 : error("all arguments should have same number of values");
259 : }
260 20 : if( getPntrToArgument(i)->getRank()==1 ) {
261 14 : argnames.push_back( getPntrToArgument(i)->getName() );
262 6 : } else if( getPntrToArgument(i)->getRank()==2 ) {
263 6 : (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames );
264 : }
265 : }
266 15 : }
267 :
268 25 : void DumpPDB::update() {
269 25 : OFile opdbf;
270 25 : opdbf.link(*this);
271 25 : opdbf.setBackupString("analysis");
272 25 : opdbf.open( file );
273 25 : std::size_t psign=fmt.find("%");
274 25 : plumed_assert( psign!=std::string::npos );
275 25 : std::string descr2="%s=%-" + fmt.substr(psign+1) + " ";
276 :
277 : // the positions of the atoms (if any) are the last argument
278 25 : unsigned nargs = getNumberOfArguments(), atomarg=0;
279 25 : if( pdb_atom_indices.size()>0 ) {
280 15 : atomarg = nargs - 1;
281 : nargs = nargs - 1;
282 : }
283 : unsigned totargs = 0;
284 25 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
285 68 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
286 43 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
287 0 : error("all arguments should have same number of values");
288 : }
289 : }
290 : // the names of the arguments (but not of the positions) are rebuilt if the shapes changed
291 53 : for(unsigned i=0; i<nargs; ++i) {
292 28 : if( getPntrToArgument(i)->getRank()==1 ) {
293 19 : totargs += 1;
294 9 : } else if( getPntrToArgument(i)->getRank()==2 ) {
295 9 : totargs += getPntrToArgument(i)->getShape()[1];
296 : }
297 : }
298 25 : if( totargs!=argnames.size() ) {
299 0 : buildArgnames();
300 : }
301 :
302 25 : if( description.length()>0 ) {
303 17 : opdbf.printf("# %s AT STEP %lld TIME %f \n", description.c_str(), getStep(), getTime() );
304 : }
305 508 : for(unsigned i=0; i<nvals; ++i) {
306 : unsigned n=0;
307 1323 : for(unsigned j=0; j<nargs; j++) {
308 : Value* thisarg=getPntrToArgument(j);
309 840 : opdbf.printf("REMARK ");
310 840 : if( thisarg->getRank()==1 ) {
311 405 : opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i) );
312 405 : n++;
313 435 : } else if( thisarg->getRank()==2 ) {
314 435 : unsigned ncols = thisarg->getShape()[1];
315 1297 : for(unsigned k=0; k<ncols; ++k) {
316 862 : opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i*ncols+k) );
317 862 : n++;
318 : }
319 : }
320 840 : opdbf.printf("\n");
321 : }
322 483 : if( pdb_atom_indices.size()>0 ) {
323 192 : unsigned npos = pdb_atom_indices.size();
324 192 : Vector pos;
325 3037 : for(unsigned k=0; k<npos; ++k) {
326 2845 : pos[0]=getPntrToArgument(atomarg)->get(npos*(3*i+0) + k);
327 2845 : pos[1]=getPntrToArgument(atomarg)->get(npos*(3*i+1) + k);
328 2845 : pos[2]=getPntrToArgument(atomarg)->get(npos*(3*i+2) + k);
329 2845 : printAtom( opdbf, pdb_atom_indices[k].serial(), pdb_resid_indices[k], pos, occupancy[k], beta[k] );
330 : }
331 : }
332 483 : opdbf.printf("END\n");
333 : }
334 25 : opdbf.close();
335 25 : }
336 :
337 : }
338 : }
|
