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 7 : void calculate() override {}
132 7 : void apply() override {}
133 : void update() override ;
134 : };
135 :
136 : PLUMED_REGISTER_ACTION(DumpPDB,"DUMPPDB")
137 :
138 21 : void DumpPDB::registerKeywords( Keywords& keys ) {
139 21 : Action::registerKeywords( keys );
140 21 : ActionPilot::registerKeywords( keys );
141 21 : ActionWithArguments::registerKeywords( keys );
142 42 : keys.addInputKeyword("optional","ARG","vector/matrix","the values that are being output in the PDB file");
143 21 : keys.add("optional","ATOMS","value containing positions of atoms that should be output");
144 21 : keys.add("compulsory","STRIDE","0","the frequency with which the atoms should be output");
145 21 : keys.add("hidden","FMT","the format that should be used to output real numbers in the title");
146 21 : keys.add("compulsory","FILE","the name of the file on which to output these quantities");
147 21 : keys.add("compulsory","OCCUPANCY","1.0","vector of values to output in the occupancy column of the pdb file");
148 21 : keys.add("compulsory","BETA","1.0","vector of values to output in the beta column of the pdb file");
149 21 : keys.add("optional","DESCRIPTION","the title to use for your PDB output");
150 21 : keys.add("optional","ATOM_INDICES","the indices of the atoms in your PDB output");
151 21 : keys.add("optional","RESIDUE_INDICES","the indices of the residues in your PDB output");
152 21 : keys.add("optional","ARG_NAMES","the names of the arguments that are being output");
153 21 : }
154 :
155 15 : DumpPDB::DumpPDB(const ActionOptions&ao):
156 : Action(ao),
157 : ActionWithArguments(ao),
158 15 : ActionPilot(ao) {
159 30 : parse("FILE",file);
160 15 : if(file.length()==0) {
161 0 : error("name out output file was not specified");
162 : }
163 : std::vector<std::string> atoms;
164 30 : parseVector("ATOMS",atoms);
165 15 : if( atoms.size()>0 ) {
166 : std::vector<Value*> atom_args;
167 8 : interpretArgumentList( atoms, plumed.getActionSet(), this, atom_args );
168 8 : if( atom_args.size()!=1 ) {
169 0 : error("only one action should appear in input to atom args");
170 : }
171 8 : std::vector<Value*> args( getArguments() );
172 8 : args.push_back( atom_args[0] );
173 8 : requestArguments( args );
174 : std::vector<std::string> indices;
175 16 : parseVector("ATOM_INDICES",indices);
176 : std::vector<Value*> xvals,yvals,zvals,masv,chargev;
177 8 : ActionAtomistic::getAtomValuesFromPlumedObject(plumed,xvals,yvals,zvals,masv,chargev);
178 8 : ActionAtomistic::interpretAtomList( indices, xvals, this, pdb_atom_indices );
179 8 : log.printf(" printing atoms : ");
180 153 : for(unsigned i=0; i<indices.size(); ++i) {
181 145 : log.printf("%d ", pdb_atom_indices[i].serial() );
182 : }
183 8 : log.printf("\n");
184 8 : 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 8 : beta.resize( atom_args[0]->getShape()[1]/3 );
188 8 : occupancy.resize( atom_args[0]->getShape()[1]/3 );
189 8 : parseVector("OCCUPANCY", occupancy );
190 8 : parseVector("BETA", beta );
191 16 : parseVector("RESIDUE_INDICES",pdb_resid_indices);
192 8 : 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 3 : } 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 8 : }
201 15 : log.printf(" printing configurations to PDB file to file named %s \n", file.c_str() );
202 30 : parseVector("ARG_NAMES",argnames);
203 15 : if( argnames.size()==0 ) {
204 14 : buildArgnames();
205 : }
206 : fmt="%f";
207 15 : parse("FMT",fmt);
208 15 : fmt=" "+fmt;
209 15 : if( getStride()==0 ) {
210 8 : setStride(0);
211 8 : log.printf(" printing pdb at end of calculation \n");
212 : }
213 :
214 30 : parse("DESCRIPTION",description);
215 15 : if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->hasDerivatives() ) {
216 0 : error("argument for printing of PDB should be vector or matrix");
217 : }
218 15 : unsigned nv=getPntrToArgument(0)->getShape()[0];
219 44 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
220 29 : if( getPntrToArgument(i)->getRank()==0 || getPntrToArgument(i)->hasDerivatives() ) {
221 0 : error("argument for printing of PDB should be vector or matrix");
222 : }
223 29 : if( getPntrToArgument(i)->getShape()[0]!=nv ) {
224 0 : error("for printing to pdb files all arguments must have same number of values");
225 : }
226 : }
227 15 : }
228 :
229 2143 : void DumpPDB::printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const {
230 2143 : if( rnum>999 ) {
231 0 : plumed_merror("atom number too large to be used as residue number");
232 : }
233 : std::array<char,6> at;
234 2143 : const char* msg = h36::hy36encode(5,anum,&at[0]);
235 2143 : plumed_assert(msg==nullptr) << msg;
236 2143 : at[5]=0;
237 2143 : double lunits = plumed.getUnits().getLength()/0.1;
238 2143 : opdbf.printf("ATOM %s X RES %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
239 : &at[0], rnum, lunits*pos[0], lunits*pos[1], lunits*pos[2], m, q );
240 2143 : }
241 :
242 23 : void DumpPDB::buildArgnames() {
243 23 : unsigned nargs = getNumberOfArguments();
244 23 : if( pdb_atom_indices.size()>0 ) {
245 17 : nargs = nargs - 1;
246 : }
247 : if( nargs<0 ) {
248 : return;
249 : }
250 :
251 23 : argnames.resize(0);
252 23 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
253 : if( getPntrToArgument(0)->getRank()==2 ) {
254 : nvals = getPntrToArgument(0)->getShape()[0];
255 : }
256 48 : for(unsigned i=0; i<nargs; ++i) {
257 25 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
258 0 : error("all arguments should have same number of values");
259 : }
260 25 : if( getPntrToArgument(i)->getRank()==1 ) {
261 17 : argnames.push_back( getPntrToArgument(i)->getName() );
262 8 : } else if( getPntrToArgument(i)->getRank()==2 ) {
263 8 : (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames );
264 : }
265 : }
266 : }
267 :
268 19 : void DumpPDB::update() {
269 19 : OFile opdbf;
270 19 : opdbf.link(*this);
271 19 : opdbf.setBackupString("analysis");
272 19 : opdbf.open( file );
273 19 : std::size_t psign=fmt.find("%");
274 19 : plumed_assert( psign!=std::string::npos );
275 38 : std::string descr2="%s=%-" + fmt.substr(psign+1) + " ";
276 :
277 : unsigned totargs = 0;
278 19 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
279 56 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
280 37 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
281 0 : error("all arguments should have same number of values");
282 : }
283 37 : if( getPntrToArgument(i)->getRank()==1 ) {
284 19 : totargs += 1;
285 18 : } else if( getPntrToArgument(i)->getRank()==2 ) {
286 18 : totargs += getPntrToArgument(i)->getShape()[1];
287 : }
288 : }
289 19 : if( totargs!=argnames.size() ) {
290 9 : buildArgnames();
291 : }
292 :
293 19 : if( description.length()>0 ) {
294 11 : opdbf.printf("# %s AT STEP %lld TIME %f \n", description.c_str(), getStep(), getTime() );
295 : }
296 19 : unsigned nargs = getNumberOfArguments(), atomarg=0;
297 19 : if( pdb_atom_indices.size()>0 ) {
298 9 : if( nargs>1 ) {
299 3 : atomarg = nargs - 1;
300 : nargs = nargs-1;
301 : } else {
302 : nargs = 0;
303 : }
304 : }
305 448 : for(unsigned i=0; i<nvals; ++i) {
306 : unsigned n=0;
307 1269 : 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 429 : if( pdb_atom_indices.size()>0 ) {
323 138 : unsigned npos = pdb_atom_indices.size();
324 : Vector pos;
325 2281 : for(unsigned k=0; k<npos; ++k) {
326 2143 : pos[0]=getPntrToArgument(atomarg)->get(npos*(3*i+0) + k);
327 2143 : pos[1]=getPntrToArgument(atomarg)->get(npos*(3*i+1) + k);
328 2143 : pos[2]=getPntrToArgument(atomarg)->get(npos*(3*i+2) + k);
329 2143 : printAtom( opdbf, pdb_atom_indices[k].serial(), pdb_resid_indices[k], pos, occupancy[k], beta[k] );
330 : }
331 : }
332 429 : opdbf.printf("END\n");
333 : }
334 19 : opdbf.close();
335 19 : }
336 :
337 : }
338 : }
|
