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 : Notice, furthermore, that you can output the positions of atoms along with the argument values that correspond to each
52 : set of atomic positions as follows:
53 :
54 : ```plumed
55 : # Calculate the distance between atoms 1 and 2
56 : d: DISTANCE ATOMS=1,2
57 : # Collect the distance between atoms 1 and 2
58 : cd: COLLECT ARG=d
59 : # Calculate the angle involving atoms 1, 2 and 3
60 : a: ANGLE ATOMS=1,2,3
61 : # Collect the angle involving atoms 1, 2 and 3
62 : ca: COLLECT ARG=a
63 : # Now collect the positions
64 : ff: COLLECT_FRAMES ATOMS=1-22 STRIDE=1
65 : # Output a PDB file that contains the collected atomic positions
66 : # and the collected distances and angles
67 : DUMPPDB ATOM_INDICES=1-22 ATOMS=ff_data ARG=cd,ca FILE=pos.pdb
68 : ```
69 :
70 : The outputted pdb file has the format described in the documentation for the [PDB2CONSTANT](PDB2CONSTANT.md) action.
71 :
72 : */
73 : //+ENDPLUMEDOC
74 :
75 : class DumpPDB :
76 : public ActionWithArguments,
77 : public ActionPilot {
78 : std::string fmt;
79 : std::string file;
80 : std::string description;
81 : std::vector<unsigned> pdb_resid_indices;
82 : std::vector<double> beta, occupancy;
83 : std::vector<std::string> argnames;
84 : std::vector<AtomNumber> pdb_atom_indices;
85 : void buildArgnames();
86 : void printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const ;
87 : public:
88 : static void registerKeywords( Keywords& keys );
89 : explicit DumpPDB(const ActionOptions&);
90 7 : void calculate() override {}
91 7 : void apply() override {}
92 : void update() override ;
93 : };
94 :
95 : PLUMED_REGISTER_ACTION(DumpPDB,"DUMPPDB")
96 :
97 21 : void DumpPDB::registerKeywords( Keywords& keys ) {
98 21 : Action::registerKeywords( keys );
99 21 : ActionPilot::registerKeywords( keys );
100 21 : ActionWithArguments::registerKeywords( keys );
101 42 : keys.addInputKeyword("optional","ARG","vector/matrix","the values that are being output in the PDB file");
102 21 : keys.add("optional","ATOMS","value containing positions of atoms that should be output");
103 21 : keys.add("compulsory","STRIDE","0","the frequency with which the atoms should be output");
104 21 : keys.add("compulsory","FILE","the name of the file on which to output these quantities");
105 21 : keys.add("compulsory","FMT","%f","the format that should be used to output real numbers");
106 21 : keys.add("compulsory","OCCUPANCY","1.0","vector of values to output in the occupancy column of the pdb file");
107 21 : keys.add("compulsory","BETA","1.0","vector of values to output in the beta column of the pdb file");
108 21 : keys.add("optional","DESCRIPTION","the title to use for your PDB output");
109 21 : keys.add("optional","ATOM_INDICES","the indices of the atoms in your PDB output");
110 21 : keys.add("optional","RESIDUE_INDICES","the indices of the residues in your PDB output");
111 21 : keys.add("optional","ARG_NAMES","the names of the arguments that are being output");
112 21 : }
113 :
114 15 : DumpPDB::DumpPDB(const ActionOptions&ao):
115 : Action(ao),
116 : ActionWithArguments(ao),
117 15 : ActionPilot(ao) {
118 30 : parse("FILE",file);
119 15 : if(file.length()==0) {
120 0 : error("name out output file was not specified");
121 : }
122 : std::vector<std::string> atoms;
123 30 : parseVector("ATOMS",atoms);
124 15 : if( atoms.size()>0 ) {
125 : std::vector<Value*> atom_args;
126 8 : interpretArgumentList( atoms, plumed.getActionSet(), this, atom_args );
127 8 : if( atom_args.size()!=1 ) {
128 0 : error("only one action should appear in input to atom args");
129 : }
130 8 : std::vector<Value*> args( getArguments() );
131 8 : args.push_back( atom_args[0] );
132 8 : requestArguments( args );
133 : std::vector<std::string> indices;
134 16 : parseVector("ATOM_INDICES",indices);
135 : std::vector<Value*> xvals,yvals,zvals,masv,chargev;
136 8 : ActionAtomistic::getAtomValuesFromPlumedObject(plumed,xvals,yvals,zvals,masv,chargev);
137 8 : ActionAtomistic::interpretAtomList( indices, xvals, this, pdb_atom_indices );
138 8 : log.printf(" printing atoms : ");
139 153 : for(unsigned i=0; i<indices.size(); ++i) {
140 145 : log.printf("%d ", pdb_atom_indices[i].serial() );
141 : }
142 8 : log.printf("\n");
143 8 : if( pdb_atom_indices.size()!=atom_args[0]->getShape()[1]/3 ) {
144 0 : error("mismatch between size of matrix containing positions and vector of atom indices");
145 : }
146 8 : beta.resize( atom_args[0]->getShape()[1]/3 );
147 8 : occupancy.resize( atom_args[0]->getShape()[1]/3 );
148 8 : parseVector("OCCUPANCY", occupancy );
149 8 : parseVector("BETA", beta );
150 16 : parseVector("RESIDUE_INDICES",pdb_resid_indices);
151 8 : if( pdb_resid_indices.size()==0 ) {
152 5 : pdb_resid_indices.resize( pdb_atom_indices.size() );
153 111 : for(unsigned i=0; i<pdb_resid_indices.size(); ++i) {
154 106 : pdb_resid_indices[i] = pdb_atom_indices[i].serial();
155 : }
156 3 : } else if( pdb_resid_indices.size()!=pdb_atom_indices.size() ) {
157 0 : error("mismatch between number of atom indices provided and number of residue indices provided");
158 : }
159 8 : }
160 15 : log.printf(" printing configurations to PDB file to file named %s \n", file.c_str() );
161 30 : parseVector("ARG_NAMES",argnames);
162 15 : if( argnames.size()==0 ) {
163 14 : buildArgnames();
164 : }
165 15 : parse("FMT",fmt);
166 15 : fmt=" "+fmt;
167 15 : if( getStride()==0 ) {
168 8 : setStride(0);
169 8 : log.printf(" printing pdb at end of calculation \n");
170 : }
171 :
172 30 : parse("DESCRIPTION",description);
173 15 : if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->hasDerivatives() ) {
174 0 : error("argument for printing of PDB should be vector or matrix");
175 : }
176 15 : unsigned nv=getPntrToArgument(0)->getShape()[0];
177 44 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
178 29 : if( getPntrToArgument(i)->getRank()==0 || getPntrToArgument(i)->hasDerivatives() ) {
179 0 : error("argument for printing of PDB should be vector or matrix");
180 : }
181 29 : if( getPntrToArgument(i)->getShape()[0]!=nv ) {
182 0 : error("for printing to pdb files all arguments must have same number of values");
183 : }
184 : }
185 15 : }
186 :
187 2143 : void DumpPDB::printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const {
188 2143 : if( rnum>999 ) {
189 0 : plumed_merror("atom number too large to be used as residue number");
190 : }
191 : std::array<char,6> at;
192 2143 : const char* msg = h36::hy36encode(5,anum,&at[0]);
193 2143 : plumed_assert(msg==nullptr) << msg;
194 2143 : at[5]=0;
195 2143 : double lunits = plumed.getUnits().getLength()/0.1;
196 2143 : opdbf.printf("ATOM %s X RES %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
197 2143 : &at[0], rnum, lunits*pos[0], lunits*pos[1], lunits*pos[2], m, q );
198 2143 : }
199 :
200 23 : void DumpPDB::buildArgnames() {
201 23 : unsigned nargs = getNumberOfArguments();
202 23 : if( pdb_atom_indices.size()>0 ) {
203 17 : nargs = nargs - 1;
204 : }
205 : if( nargs<0 ) {
206 : return;
207 : }
208 :
209 23 : argnames.resize(0);
210 23 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
211 : if( getPntrToArgument(0)->getRank()==2 ) {
212 : nvals = getPntrToArgument(0)->getShape()[0];
213 : }
214 48 : for(unsigned i=0; i<nargs; ++i) {
215 25 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
216 0 : error("all arguments should have same number of values");
217 : }
218 25 : if( getPntrToArgument(i)->getRank()==1 ) {
219 17 : argnames.push_back( getPntrToArgument(i)->getName() );
220 8 : } else if( getPntrToArgument(i)->getRank()==2 ) {
221 8 : (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames );
222 : }
223 25 : getPntrToArgument(i)->buildDataStore();
224 : }
225 : }
226 :
227 19 : void DumpPDB::update() {
228 19 : OFile opdbf;
229 19 : opdbf.link(*this);
230 19 : opdbf.setBackupString("analysis");
231 19 : opdbf.open( file );
232 19 : std::size_t psign=fmt.find("%");
233 19 : plumed_assert( psign!=std::string::npos );
234 38 : std::string descr2="%s=%-" + fmt.substr(psign+1) + " ";
235 :
236 : unsigned totargs = 0;
237 19 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
238 56 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
239 37 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
240 0 : error("all arguments should have same number of values");
241 : }
242 37 : if( getPntrToArgument(i)->getRank()==1 ) {
243 19 : totargs += 1;
244 18 : } else if( getPntrToArgument(i)->getRank()==2 ) {
245 18 : totargs += getPntrToArgument(i)->getShape()[1];
246 : }
247 : }
248 19 : if( totargs!=argnames.size() ) {
249 9 : buildArgnames();
250 : }
251 :
252 19 : if( description.length()>0 ) {
253 11 : opdbf.printf("# %s AT STEP %lld TIME %f \n", description.c_str(), getStep(), getTime() );
254 : }
255 19 : unsigned nargs = getNumberOfArguments(), atomarg=0;
256 19 : if( pdb_atom_indices.size()>0 ) {
257 9 : if( nargs>1 ) {
258 3 : atomarg = nargs - 1;
259 : nargs = nargs-1;
260 : } else {
261 : nargs = 0;
262 : }
263 : }
264 448 : for(unsigned i=0; i<nvals; ++i) {
265 : unsigned n=0;
266 1269 : for(unsigned j=0; j<nargs; j++) {
267 : Value* thisarg=getPntrToArgument(j);
268 840 : opdbf.printf("REMARK ");
269 840 : if( thisarg->getRank()==1 ) {
270 405 : opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i) );
271 405 : n++;
272 435 : } else if( thisarg->getRank()==2 ) {
273 435 : unsigned ncols = thisarg->getShape()[1];
274 1297 : for(unsigned k=0; k<ncols; ++k) {
275 862 : opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i*ncols+k) );
276 862 : n++;
277 : }
278 : }
279 840 : opdbf.printf("\n");
280 : }
281 429 : if( pdb_atom_indices.size()>0 ) {
282 138 : unsigned npos = pdb_atom_indices.size();
283 138 : Vector pos;
284 2281 : for(unsigned k=0; k<npos; ++k) {
285 2143 : pos[0]=getPntrToArgument(atomarg)->get(npos*(3*i+0) + k);
286 2143 : pos[1]=getPntrToArgument(atomarg)->get(npos*(3*i+1) + k);
287 2143 : pos[2]=getPntrToArgument(atomarg)->get(npos*(3*i+2) + k);
288 2143 : printAtom( opdbf, pdb_atom_indices[k].serial(), pdb_resid_indices[k], pos, occupancy[k], beta[k] );
289 : }
290 : }
291 429 : opdbf.printf("END\n");
292 : }
293 19 : opdbf.close();
294 19 : }
295 :
296 : }
297 : }
|
