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/ActionAtomistic.h"
23 : #include "core/ActionWithArguments.h"
24 : #include "core/ActionPilot.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Pbc.h"
27 : #include "tools/File.h"
28 : #include "core/PlumedMain.h"
29 : #include "tools/Units.h"
30 : #include "tools/CheckInRange.h"
31 : #include <cstdio>
32 : #include <memory>
33 : #include "core/GenericMolInfo.h"
34 : #include "core/ActionSet.h"
35 : #include "xdrfile/xdrfile_xtc.h"
36 : #include "xdrfile/xdrfile_trr.h"
37 :
38 :
39 : namespace PLMD {
40 : namespace generic {
41 :
42 : //+PLUMEDOC PRINTANALYSIS DUMPATOMS
43 : /*
44 : Dump selected atoms on a file.
45 :
46 : This command can be used to output the positions of a particular set of atoms.
47 : The atoms required are output in a xyz or gro formatted file.
48 : If PLUMED has been compiled with xdrfile support, then also xtc and trr files can be written.
49 : To this aim one should install xdrfile library (http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library).
50 : If the xdrfile library is installed properly the PLUMED configure script should be able to
51 : detect it and enable it.
52 : The type of file is automatically detected from the file extension, but can be also
53 : enforced with TYPE.
54 : Importantly, if your
55 : input file contains actions that edit the atoms position (e.g. \ref WHOLEMOLECULES)
56 : and the DUMPATOMS command appears after this instruction, then the edited
57 : atom positions are output.
58 : You can control the buffering of output using the \ref FLUSH keyword on a separate line.
59 :
60 : Units of the printed file can be controlled with the UNITS keyword. By default PLUMED units as
61 : controlled in the \ref UNITS command are used, but one can override it e.g. with UNITS=A.
62 : Notice that gro/xtc/trr files can only contain coordinates in nm.
63 :
64 : \par Examples
65 :
66 : The following input instructs plumed to print out the positions of atoms
67 : 1-10 together with the position of the center of mass of atoms 11-20 every
68 : 10 steps to a file called file.xyz.
69 : \plumedfile
70 : COM ATOMS=11-20 LABEL=c1
71 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1
72 : \endplumedfile
73 : Notice that the coordinates in the xyz file will be expressed in nm, since these
74 : are the defaults units in PLUMED. If you want the xyz file to be expressed in A, you should use the
75 : following input
76 : \plumedfile
77 : COM ATOMS=11-20 LABEL=c1
78 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1 UNITS=A
79 : \endplumedfile
80 : As an alternative, you might want to set all the length used by PLUMED to Angstrom using the \ref UNITS
81 : action. However, this latter choice will affect all your input and output.
82 :
83 : The following input is very similar but dumps a .gro (gromacs) file,
84 : which also contains atom and residue names.
85 : \plumedfile
86 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
87 : # this is required to have proper atom names:
88 : MOLINFO STRUCTURE=reference.pdb
89 : # if omitted, atoms will have "X" name...
90 :
91 : COM ATOMS=11-20 LABEL=c1
92 : DUMPATOMS STRIDE=10 FILE=file.gro ATOMS=1-10,c1
93 : # notice that last atom is a virtual one and will not have
94 : # a correct name in the resulting gro file
95 : \endplumedfile
96 :
97 : The `file.gro` will contain coordinates expressed in nm, since this is the convention for gro files.
98 :
99 : In case you have compiled PLUMED with `xdrfile` library, you might even write xtc or trr files as follows
100 : \plumedfile
101 : COM ATOMS=11-20 LABEL=c1
102 : DUMPATOMS STRIDE=10 FILE=file.xtc ATOMS=1-10,c1
103 : \endplumedfile
104 : Notice that xtc files are significantly smaller than gro and xyz files.
105 :
106 : Finally, consider that gro and xtc file store coordinates with limited precision set by the
107 : `PRECISION` keyword. Default value is 3, which means "3 digits after dot" in nm (1/1000 of a nm).
108 : The following will write a larger xtc file with high resolution coordinates:
109 : \plumedfile
110 : COM ATOMS=11-20 LABEL=c1
111 : DUMPATOMS STRIDE=10 FILE=file.xtc ATOMS=1-10,c1 PRECISION=7
112 : \endplumedfile
113 :
114 :
115 :
116 : */
117 : //+ENDPLUMEDOC
118 :
119 : class DumpAtoms:
120 : public ActionAtomistic,
121 : public ActionWithArguments,
122 : public ActionPilot {
123 : OFile of;
124 : double lenunit;
125 : int iprecision;
126 : CheckInRange bounds;
127 : std::vector<std::string> names;
128 : std::vector<unsigned> residueNumbers;
129 : std::vector<std::string> residueNames;
130 : std::string type;
131 : std::string fmt_gro_pos;
132 : std::string fmt_gro_box;
133 : std::string fmt_xyz;
134 : xdrfile::XDRFILE* xd;
135 : public:
136 : explicit DumpAtoms(const ActionOptions&);
137 : ~DumpAtoms();
138 : static void registerKeywords( Keywords& keys );
139 : void calculateNumericalDerivatives( ActionWithValue* a=NULL ) override;
140 194 : bool actionHasForces() override {
141 194 : return false;
142 : }
143 : void lockRequests() override;
144 : void unlockRequests() override;
145 9095 : void calculate() override {}
146 9095 : void apply() override {}
147 : void update() override ;
148 : };
149 :
150 : PLUMED_REGISTER_ACTION(DumpAtoms,"DUMPATOMS")
151 :
152 226 : void DumpAtoms::registerKeywords( Keywords& keys ) {
153 226 : Action::registerKeywords( keys );
154 226 : ActionPilot::registerKeywords( keys );
155 226 : ActionAtomistic::registerKeywords( keys );
156 226 : ActionWithArguments::registerKeywords( keys );
157 226 : keys.use("ARG");
158 452 : keys.add("compulsory","STRIDE","1","the frequency with which the atoms should be output");
159 452 : keys.add("atoms", "ATOMS", "the atom indices whose positions you would like to print out");
160 452 : keys.add("compulsory", "FILE", "file on which to output coordinates; extension is automatically detected");
161 452 : keys.add("compulsory", "UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
162 452 : keys.add("optional", "PRECISION","The number of digits in trajectory file");
163 452 : keys.add("optional", "TYPE","file type, either xyz, gro, xtc, or trr, can override an automatically detected file extension");
164 452 : keys.add("optional","LESS_THAN_OR_EQUAL","when printing with arguments that are vectors only print components of vectors have a value less than or equal to this value");
165 452 : keys.add("optional","GREATER_THAN_OR_EQUAL","when printing with arguments that are vectors only print components of vectors have a value greater than or equal to this value");
166 226 : keys.use("RESTART");
167 226 : keys.use("UPDATE_FROM");
168 226 : keys.use("UPDATE_UNTIL");
169 226 : }
170 :
171 222 : DumpAtoms::DumpAtoms(const ActionOptions&ao):
172 : Action(ao),
173 : ActionAtomistic(ao),
174 : ActionWithArguments(ao),
175 : ActionPilot(ao),
176 222 : iprecision(3) {
177 : std::vector<AtomNumber> atoms;
178 : std::string file;
179 444 : parse("FILE",file);
180 222 : if(file.length()==0) {
181 0 : error("name out output file was not specified");
182 : }
183 222 : type=Tools::extension(file);
184 222 : log<<" file name "<<file<<"\n";
185 461 : if(type=="gro" || type=="xyz" || type=="xtc" || type=="trr") {
186 197 : log<<" file extension indicates a "<<type<<" file\n";
187 : } else {
188 25 : log<<" file extension not detected, assuming xyz\n";
189 : type="xyz";
190 : }
191 : std::string ntype;
192 444 : parse("TYPE",ntype);
193 222 : if(ntype.length()>0) {
194 5 : if(ntype!="xyz" && ntype!="gro" && ntype!="xtc" && ntype!="trr"
195 : ) {
196 0 : error("TYPE cannot be understood");
197 : }
198 2 : log<<" file type enforced to be "<<ntype<<"\n";
199 : type=ntype;
200 : }
201 :
202 : fmt_gro_pos="%8.3f";
203 : fmt_gro_box="%12.7f";
204 : fmt_xyz="%f";
205 :
206 : std::string precision;
207 444 : parse("PRECISION",precision);
208 222 : if(precision.length()>0) {
209 136 : Tools::convert(precision,iprecision);
210 136 : log<<" with precision "<<iprecision<<"\n";
211 : std::string a,b;
212 136 : Tools::convert(iprecision+5,a);
213 136 : Tools::convert(iprecision,b);
214 272 : fmt_gro_pos="%"+a+"."+b+"f";
215 : fmt_gro_box=fmt_gro_pos;
216 : fmt_xyz=fmt_gro_box;
217 : }
218 :
219 444 : parseAtomList("ATOMS",atoms);
220 :
221 : std::string unitname;
222 444 : parse("UNITS",unitname);
223 222 : if(unitname!="PLUMED") {
224 17 : Units myunit;
225 17 : myunit.setLength(unitname);
226 32 : if(myunit.getLength()!=1.0 && type=="gro") {
227 0 : error("gro files should be in nm");
228 : }
229 32 : if(myunit.getLength()!=1.0 && type=="xtc") {
230 0 : error("xtc files should be in nm");
231 : }
232 32 : if(myunit.getLength()!=1.0 && type=="trr") {
233 0 : error("trr files should be in nm");
234 : }
235 17 : lenunit=getUnits().getLength()/myunit.getLength();
236 534 : } else if(type=="gro" || type=="xtc" || type=="trr") {
237 56 : lenunit=getUnits().getLength();
238 : } else {
239 149 : lenunit=1.0;
240 : }
241 :
242 222 : of.link(*this);
243 222 : of.open(file);
244 : std::string path=of.getPath();
245 222 : log<<" Writing on file "<<path<<"\n";
246 : std::string mode=of.getMode();
247 222 : if(type=="xtc") {
248 6 : of.close();
249 6 : xd=xdrfile::xdrfile_open(path.c_str(),mode.c_str());
250 216 : } else if(type=="trr") {
251 4 : of.close();
252 4 : xd=xdrfile::xdrfile_open(path.c_str(),mode.c_str());
253 : }
254 222 : log.printf(" printing the following atoms in %s :", unitname.c_str() );
255 30874 : for(unsigned i=0; i<atoms.size(); ++i) {
256 30652 : log.printf(" %d",atoms[i].serial() );
257 : }
258 222 : log.printf("\n");
259 :
260 222 : if( getNumberOfArguments()>0 ) {
261 39 : if( type!="xyz" ) {
262 0 : error("can only print atomic properties when outputting xyz files");
263 : }
264 :
265 : std::vector<std::string> argnames;
266 119 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
267 80 : if( getPntrToArgument(i)->getRank()!=1 || getPntrToArgument(i)->hasDerivatives() ) {
268 0 : error("arguments for xyz output should be vectors");
269 : }
270 80 : if( getPntrToArgument(i)->getNumberOfValues()!=atoms.size() ) {
271 0 : error("number of elements in vector " + getPntrToArgument(i)->getName() + " is not equal to number of atoms output");
272 : }
273 80 : getPntrToArgument(i)->buildDataStore(true);
274 80 : argnames.push_back( getPntrToArgument(i)->getName() );
275 : }
276 : std::vector<std::string> str_upper, str_lower;
277 : std::string errors;
278 39 : parseVector("LESS_THAN_OR_EQUAL",str_upper);
279 78 : parseVector("GREATER_THAN_OR_EQUAL",str_lower);
280 39 : if( !bounds.setBounds( getNumberOfArguments(), str_lower, str_upper, errors ) ) {
281 0 : error( errors );
282 : }
283 39 : if( bounds.wereSet() ) {
284 26 : log.printf(" %s \n", bounds.report( argnames ).c_str() );
285 : }
286 39 : }
287 :
288 222 : requestAtoms(atoms, false);
289 222 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
290 222 : if( moldat ) {
291 56 : log<<" MOLINFO DATA found with label " <<moldat->getLabel()<<", using proper atom names\n";
292 56 : names.resize(atoms.size());
293 5835 : for(unsigned i=0; i<atoms.size(); i++)
294 5779 : if(atoms[i].index()<moldat->getPDBsize()) {
295 11554 : names[i]=moldat->getAtomName(atoms[i]);
296 : }
297 56 : residueNumbers.resize(atoms.size());
298 5835 : for(unsigned i=0; i<residueNumbers.size(); ++i)
299 5779 : if(atoms[i].index()<moldat->getPDBsize()) {
300 5777 : residueNumbers[i]=moldat->getResidueNumber(atoms[i]);
301 : }
302 56 : residueNames.resize(atoms.size());
303 5835 : for(unsigned i=0; i<residueNames.size(); ++i)
304 5779 : if(atoms[i].index()<moldat->getPDBsize()) {
305 11554 : residueNames[i]=moldat->getResidueName(atoms[i]);
306 : }
307 : }
308 222 : }
309 :
310 0 : void DumpAtoms::calculateNumericalDerivatives( ActionWithValue* a ) {
311 0 : plumed_merror("this should never be called");
312 : }
313 :
314 9095 : void DumpAtoms::lockRequests() {
315 : ActionWithArguments::lockRequests();
316 : ActionAtomistic::lockRequests();
317 9095 : }
318 :
319 9095 : void DumpAtoms::unlockRequests() {
320 : ActionWithArguments::unlockRequests();
321 : ActionAtomistic::unlockRequests();
322 9095 : }
323 :
324 9000 : void DumpAtoms::update() {
325 9000 : if(type=="xyz") {
326 : unsigned nat=0;
327 8414 : std::vector<double> args( getNumberOfArguments() );
328 74007 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
329 95495 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
330 29902 : args[j] = getPntrToArgument(j)->get(i);
331 : }
332 65593 : if( bounds.check( args ) ) {
333 61787 : nat++;
334 : }
335 : }
336 8414 : of.printf("%d\n",nat);
337 8414 : const Tensor & t(getPbc().getBox());
338 8414 : if(getPbc().isOrthorombic()) {
339 612 : of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2));
340 : } else {
341 16216 : of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),
342 8108 : lenunit*t(0,0),lenunit*t(0,1),lenunit*t(0,2),
343 8108 : lenunit*t(1,0),lenunit*t(1,1),lenunit*t(1,2),
344 8108 : lenunit*t(2,0),lenunit*t(2,1),lenunit*t(2,2)
345 : );
346 : }
347 74007 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
348 95495 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
349 29902 : args[j] = getPntrToArgument(j)->get(i);
350 : }
351 65593 : if( !bounds.check(args) ) {
352 3806 : continue;
353 : }
354 : const char* defname="X";
355 : const char* name=defname;
356 61787 : if(names.size()>0)
357 10802 : if(names[i].length()>0) {
358 : name=names[i].c_str();
359 : }
360 123574 : of.printf(("%s "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz).c_str(),name,lenunit*getPosition(i)(0),lenunit*getPosition(i)(1),lenunit*getPosition(i)(2));
361 87883 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
362 52192 : of.printf((" "+fmt_xyz).c_str(), getPntrToArgument(j)->get(i) );
363 : }
364 61787 : of.printf("\n");
365 : }
366 586 : } else if(type=="gro") {
367 466 : const Tensor & t(getPbc().getBox());
368 466 : of.printf("Made with PLUMED t=%f\n",getTime()/getUnits().getTime());
369 466 : of.printf("%d\n",getNumberOfAtoms());
370 38404 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
371 : const char* defname="X";
372 : const char* name=defname;
373 : unsigned residueNumber=0;
374 37938 : if(names.size()>0)
375 37137 : if(names[i].length()>0) {
376 : name=names[i].c_str();
377 : }
378 37938 : if(residueNumbers.size()>0) {
379 37137 : residueNumber=residueNumbers[i];
380 : }
381 37938 : std::string resname="";
382 37938 : if(residueNames.size()>0) {
383 37137 : resname=residueNames[i];
384 : }
385 75876 : of.printf(("%5u%-5s%5s%5d"+fmt_gro_pos+fmt_gro_pos+fmt_gro_pos+"\n").c_str(),
386 : residueNumber%100000,resname.c_str(),name,getAbsoluteIndex(i).serial()%100000,
387 37938 : lenunit*getPosition(i)(0),lenunit*getPosition(i)(1),lenunit*getPosition(i)(2));
388 : }
389 932 : of.printf((fmt_gro_box+" "+fmt_gro_box+" "+fmt_gro_box+" "+fmt_gro_box+" "+fmt_gro_box+" "+fmt_gro_box+" "+fmt_gro_box+" "+fmt_gro_box+" "+fmt_gro_box+"\n").c_str(),
390 466 : lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2),
391 466 : lenunit*t(0,1),lenunit*t(0,2),lenunit*t(1,0),
392 466 : lenunit*t(1,2),lenunit*t(2,0),lenunit*t(2,1));
393 168 : } else if(type=="xtc" || type=="trr") {
394 : xdrfile::matrix box;
395 120 : const Tensor & t(getPbc().getBox());
396 120 : int natoms=getNumberOfAtoms();
397 120 : int step=getStep();
398 120 : float time=getTime()/getUnits().getTime();
399 120 : float precision=Tools::fastpow(10.0,iprecision);
400 480 : for(int i=0; i<3; i++)
401 1440 : for(int j=0; j<3; j++) {
402 1080 : box[i][j]=lenunit*t(i,j);
403 : }
404 : // here we cannot use a std::vector<rvec> since it does not compile.
405 : // we thus use a std::unique_ptr<rvec[]>
406 : auto pos = Tools::make_unique<xdrfile::rvec[]>(natoms);
407 6456 : for(int i=0; i<natoms; i++)
408 25344 : for(int j=0; j<3; j++) {
409 19008 : pos[i][j]=lenunit*getPosition(i)(j);
410 : }
411 120 : if(type=="xtc") {
412 72 : write_xtc(xd,natoms,step,time,box,&pos[0],precision);
413 48 : } else if(type=="trr") {
414 48 : write_trr(xd,natoms,step,time,0.0,box,&pos[0],NULL,NULL);
415 : }
416 : } else {
417 0 : plumed_merror("unknown file type "+type);
418 : }
419 9000 : }
420 :
421 444 : DumpAtoms::~DumpAtoms() {
422 222 : if(type=="xtc") {
423 6 : xdrfile_close(xd);
424 216 : } else if(type=="trr") {
425 4 : xdrfile_close(xd);
426 : }
427 1110 : }
428 :
429 :
430 : }
431 : }
|
