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/Communicator.h"
31 : #include "tools/CheckInRange.h"
32 : #include "core/GenericMolInfo.h"
33 : #include "core/ActionSet.h"
34 : #include "xdrfile/xdrfile_xtc.h"
35 : #include "xdrfile/xdrfile_trr.h"
36 :
37 :
38 : namespace PLMD {
39 : namespace generic {
40 :
41 : //+PLUMEDOC PRINTANALYSIS DUMPATOMS
42 : /*
43 : Dump selected atoms on a file.
44 :
45 : This command can be used to output the positions of a particular set of atoms.
46 : For example, the following input instructs PLUMED to print out the positions
47 : of atoms 1-10 together with the position of the center of mass of atoms 11-20 every
48 : 10 steps to a file called file.xyz.
49 :
50 : ```plumed
51 : c1: COM ATOMS=11-20
52 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1
53 : ```
54 :
55 : By default, the coordinates in the output xyz file are in nm but you can change these units
56 : by using the `UNITS` keyword as shown below:
57 :
58 : ```plumed
59 : c1: COM ATOMS=11-20
60 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1 UNITS=A
61 : ```
62 :
63 : or by using the [UNITS](UNITS.md) action as shown below:
64 :
65 : ```plumed
66 : UNITS LENGTH=A
67 : c1: COM ATOMS=11-20
68 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1
69 : ```
70 :
71 : Notice, however, that if you use the second option here all the quantitities with units of length in your input
72 : file must be provided in Angstrom and not nm.
73 :
74 : ## DUMPATOMS and WHOLEMOLECULES
75 :
76 : The commands [WHOLEMOLECULES](WHOLEMOLECULES.md), [WRAPAROUND](WRAPAROUND.md), [FIT_TO_TEMPLATE](FIT_TO_TEMPLATE.md)
77 : and [RESET_CELL](RESET_CELL.md) all edit the global positions of the atoms. If you use an input like this one:
78 :
79 : ```plumed
80 : DUMPATOMS ATOMS=1-10 FILE=file.xyz
81 : WHOLEMOLECULES ENTITY0=1-10
82 : ```
83 :
84 : then the positions of the atoms that were passed to PLUMED by the MD code are output. However, if you use an input
85 : like this one:
86 :
87 : ```plumed
88 : WHOLEMOLECULES ENTITY0=1-10
89 : DUMPATOMS ATOMS=1-10 FILE=file.xyz
90 : ```
91 :
92 : the positions outputted by the DUMPATOMS command will have been editted by the [WHOLEMOLECULES](WHOLEMOLECULES.md) command.
93 :
94 : ## Outputting other file types
95 :
96 : The extension that is given to the file specified using the `FILE` keyword determines the output file type. Hence,
97 : the following example will output a gro file rather than an xyz file:
98 :
99 : ```plumed
100 : c1: COM ATOMS=11-20
101 : DUMPATOMS STRIDE=10 FILE=file.gro ATOMS=1-10,c1
102 : ```
103 :
104 : You can also enforce the output file type by using the `TYPE` keyword as shown below:
105 :
106 : ```plumed
107 : c1: COM ATOMS=11-20
108 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1 TYPE=gro
109 : FLUSH STRIDE=1
110 : ```
111 :
112 : Notice that DUMPATOMS command here outputs the atoms in the gro-file format even though the author of this input has used the xyz extension.
113 : Further note that by using the [FLUSH](FLUSH.md) we force PLUMED to output flush all the open files every step and not to store output
114 : data in a buffer before printing it to the output files.
115 :
116 : Outputting the atomic positions using the gro file format is particularly advantageous if you also have a [MOLINFO](MOLINFO.md) command in
117 : your input file as shown below:
118 :
119 : ```plumed
120 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
121 : # this is required to have proper atom names:
122 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
123 : # if omitted, atoms will have "X" name...
124 :
125 : c1: COM ATOMS=11-20
126 : DUMPATOMS STRIDE=10 FILE=file.gro ATOMS=1-10,c1
127 : # notice that last atom is a virtual one and will not have
128 : # a correct name in the resulting gro file
129 : ```
130 :
131 : The reason that using the gro file format is advantageous in this case is that PLUMED will also output the atom and residue names
132 : for the non-virtual atoms. PLUMED is able to do this in this case as it is able to use the information that was read in from the
133 : pdb file that was provided to the [MOLINFO](MOLINFO.md) command.
134 :
135 : If PLUMED has been compiled with xdrfile support, then PLUMED
136 : can output xtc and trr files as well as xyz and gro files. If you want to use these output types you should install the xdrfile
137 : library by following the instructions [here](http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library).
138 : If the xdrfile library is installed properly the PLUMED configure script should be able to
139 : detect it and enable it. The following example shows how you can use DUMPATOMS to output an xtc file:
140 :
141 : ```plumed
142 : c1: COM ATOMS=11-20
143 : DUMPATOMS STRIDE=10 FILE=file.xtc ATOMS=1-10,c1
144 : ```
145 :
146 : The xtc file that is output by this command will be significantly smaller than a gro and xyz file.
147 :
148 : Finally, notice that gro and xtc file store coordinates with limited precision set by the
149 : `PRECISION` keyword. The default value is 3, which means "3 digits after dot" in nm (1/1000 of a nm).
150 : The following will write a larger xtc file with high resolution coordinates:
151 :
152 : ```plumed
153 : COM ATOMS=11-20 LABEL=c1
154 : DUMPATOMS STRIDE=10 FILE=file.xtc ATOMS=1-10,c1 PRECISION=7
155 : ```
156 :
157 : ## Outputting atomic positions and vectors
158 :
159 : The atoms section of an xyz file normally contains four columns of data - a symbol that tells you the atom type
160 : and then three columns containing the atom's $x$, $y$ and $z$ coordinates. PLUMED allows you to output more columns
161 : of data in this file. For example, the following input outputs five columns of data. The first four columns of data
162 : are the usual ones you would expect in an xyz file, while the fifth contains the coordination numbers for each of the
163 : atom that have been calculated using PLUMED
164 :
165 : ```plumed
166 : # These three lines calculate the coordination numbers of 100 atoms
167 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1 NN=6 MM=12}
168 : ones: ONES SIZE=100
169 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
170 : DUMPATOMS ATOMS=1-100 ARG=cc FILE=file.xyz
171 : ```
172 :
173 : This command is used in the shortcut that recovered the old [DUMPMULTICOLVAR](DUMPMULTICOLVAR.md) command. This new version of
174 : the command is better, however, as you can output more than one vector of symmetry functions at once as is demonstrated by the following
175 : input that outputs the coordination numbers and the values that were obtained when the coordination numbers were all transformed by a
176 : switching function:
177 :
178 : ```plumed
179 : # These three lines calculate the coordination numbers of 100 atoms
180 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1 NN=6 MM=12}
181 : ones: ONES SIZE=100
182 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
183 : fcc: LESS_THAN ARG=cc SWITCH={RATIONAL R_0=4}
184 : DUMPATOMS ATOMS=1-100 ARG=cc,fcc FILE=file.xyz
185 : ```
186 :
187 : Notice that when we use an `ARG` keyword we can also use DUMPATOMS to only print those atoms whose corresponding element in the
188 : the input vector satisfies a certain criteria. For example the input file below only prints the positions (and coordination numbers) of atoms that
189 : have a coordination number that is greater than or equal to 4.
190 :
191 : ```plumed
192 : # These three lines calculate the coordination numbers of 100 atoms
193 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1 NN=6 MM=12}
194 : ones: ONES SIZE=100
195 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
196 : DUMPATOMS ATOMS=1-100 ARG=cc GREATER_THAN_OR_EQUAL=4 FILE=file.xyz
197 : ```
198 :
199 : Alternatively, the following input allows us to output those atoms that have a coordination number that is between 4 and 6:
200 :
201 : ```plumed
202 : # These three lines calculate the coordination numbers of 100 atoms
203 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1 NN=6 MM=12}
204 : ones: ONES SIZE=100
205 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
206 : DUMPATOMS ...
207 : ATOMS=1-100 ARG=cc
208 : GREATER_THAN_OR_EQUAL=4
209 : LESS_THAN_OR_EQUAL=6
210 : FILE=file.xyz
211 : ...
212 : ```
213 :
214 : Commands like these are useful if you want to print the coordinates of the atom that are in a paricular cluster that has been identified using
215 : the [DFSCLUSTERING](DFSCLUSTERING.md) command.
216 :
217 : __You can only use the ARG keyword if you are outputting an xyz file__
218 :
219 : ## PRINT and RESTART
220 :
221 : If you run a calculation with the following input:
222 :
223 : ```plumed
224 : DUMPATOMS ATOMS=1-100 FILE=file.xyz
225 : ```
226 :
227 : and a file called `file.xyz` is already present in the directory where the calculation is running, the existing file is backed up
228 : and renamed to `bck.0.file.xyz` so that new data can be output to a new file called `file.xyz`. If you would like to append to the
229 : existing file you can use the RESTART command as shown below:
230 :
231 : ```plumed
232 : DUMPATOMS ATOMS=1-100 FILE=file.xyz RESTART=YES
233 : ```
234 :
235 : You can achieve the same result by using the [RESTART](RESTART.md) action as shown below:
236 :
237 : ```plumed
238 : RESTART
239 : DUMPATOMS ATOMS=1-100 FILE=file.xyz
240 : ```
241 :
242 : However, the advantage of using the RESTART keyword is that you can apped to some files and back up others as illustrated below:
243 :
244 : ```plumed
245 : DUMPATOMS ATOMS=1-100 FILE=file1.xyz
246 : DUMPATOMS ATOMS=1-100 FILE=file2.xyz RESTART=YES
247 : ```
248 :
249 : If you use the input above the file `file1.xyz` is backed up, while new data will be appended to the file `file2.xyz`. If you use the
250 : [RESTART](RESTART.md) action instead data will be appended to both colvar files.
251 :
252 : ## Switching printing on and off
253 :
254 : You can use the UPDATE_FROM and UPDATE_UNTIL flags to make the DUMPATOMS command only output data at certain points during the trajectory.
255 : To see how this works consider the following example:
256 :
257 : ```plumed
258 : DUMPATOMS ATOMS=1-100 FILE=file.xyz UPDATE_FROM=100 UPDATE_UNTIL=500 STRIDE=1
259 : ```
260 :
261 : During the first 100 ps of a simulation with this input the atomic positions are not output to the file called file.xyz.
262 : The positions are instead only output after the first 100 ps of trajectory have elapsed. Furthermore, output of the positions stops
263 : once the trajectory is longer than 500 ps. In other words, the positions are only output during the 400 ps time interval after the first
264 : 100 ps of the simulation.
265 :
266 : */
267 : //+ENDPLUMEDOC
268 :
269 : class DumpAtoms:
270 : public ActionAtomistic,
271 : public ActionWithArguments,
272 : public ActionPilot {
273 : OFile of;
274 : double lenunit;
275 : int iprecision;
276 : CheckInRange bounds;
277 : std::vector<std::string> names;
278 : std::vector<int> residueNumbers;
279 : std::vector<std::string> residueNames;
280 : std::string type;
281 : std::string fmt_gro_pos;
282 : std::string fmt_gro_box;
283 : std::string fmt_xyz;
284 : xdrfile::XDRFILE* xd=nullptr;
285 : public:
286 : explicit DumpAtoms(const ActionOptions&);
287 : ~DumpAtoms();
288 : static void registerKeywords( Keywords& keys );
289 : void calculateNumericalDerivatives( ActionWithValue* a=NULL ) override;
290 8818 : bool actionHasForces() override {
291 8818 : return false;
292 : }
293 : void lockRequests() override;
294 : void unlockRequests() override;
295 9113 : void calculate() override {}
296 9113 : void apply() override {}
297 : void update() override ;
298 : };
299 :
300 : PLUMED_REGISTER_ACTION(DumpAtoms,"DUMPATOMS")
301 :
302 231 : void DumpAtoms::registerKeywords( Keywords& keys ) {
303 231 : Action::registerKeywords( keys );
304 231 : ActionPilot::registerKeywords( keys );
305 231 : ActionAtomistic::registerKeywords( keys );
306 231 : ActionWithArguments::registerKeywords( keys );
307 462 : keys.addInputKeyword("optional","ARG","vector","the labels of vectors that should be output in the xyz file. The number of elements in the vector should equal the number of atoms output");
308 231 : keys.add("compulsory","STRIDE","1","the frequency with which the atoms should be output");
309 231 : keys.add("atoms", "ATOMS", "the atom indices whose positions you would like to print out");
310 231 : keys.add("compulsory", "FILE", "file on which to output coordinates; extension is automatically detected");
311 231 : keys.add("compulsory", "UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
312 231 : keys.add("optional", "PRECISION","The number of digits in trajectory file");
313 231 : keys.add("optional", "TYPE","file type, either xyz, gro, xtc, or trr, can override an automatically detected file extension");
314 231 : 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");
315 231 : 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");
316 231 : keys.use("RESTART");
317 231 : keys.use("UPDATE_FROM");
318 231 : keys.use("UPDATE_UNTIL");
319 231 : }
320 :
321 227 : DumpAtoms::DumpAtoms(const ActionOptions&ao):
322 : Action(ao),
323 : ActionAtomistic(ao),
324 : ActionWithArguments(ao),
325 : ActionPilot(ao),
326 227 : iprecision(3) {
327 227 : std::vector<AtomNumber> atoms;
328 : std::string file;
329 454 : parse("FILE",file);
330 227 : if(file.length()==0) {
331 0 : error("name out output file was not specified");
332 : }
333 227 : type=Tools::extension(file);
334 227 : log<<" file name "<<file<<"\n";
335 227 : if(type=="gro" || type=="xyz" || type=="xtc" || type=="trr") {
336 202 : log<<" file extension indicates a "<<type<<" file\n";
337 : } else {
338 25 : log<<" file extension not detected, assuming xyz\n";
339 : type="xyz";
340 : }
341 : std::string ntype;
342 454 : parse("TYPE",ntype);
343 227 : if(ntype.length()>0) {
344 2 : if(ntype!="xyz" && ntype!="gro" && ntype!="xtc" && ntype!="trr"
345 : ) {
346 0 : error("TYPE cannot be understood");
347 : }
348 2 : log<<" file type enforced to be "<<ntype<<"\n";
349 : type=ntype;
350 : }
351 :
352 227 : fmt_gro_pos="%8.3f";
353 227 : fmt_gro_box="%12.7f";
354 227 : fmt_xyz="%f";
355 :
356 : std::string precision;
357 454 : parse("PRECISION",precision);
358 227 : if(precision.length()>0) {
359 141 : Tools::convert(precision,iprecision);
360 141 : log<<" with precision "<<iprecision<<"\n";
361 : std::string a,b;
362 141 : Tools::convert(iprecision+5,a);
363 141 : Tools::convert(iprecision,b);
364 564 : fmt_gro_pos="%"+a+"."+b+"f";
365 : fmt_gro_box=fmt_gro_pos;
366 : fmt_xyz=fmt_gro_box;
367 : }
368 :
369 454 : parseAtomList("ATOMS",atoms);
370 :
371 : std::string unitname;
372 227 : parse("UNITS",unitname);
373 227 : if(unitname!="PLUMED") {
374 17 : Units myunit;
375 17 : myunit.setLength(unitname);
376 17 : if(myunit.getLength()!=1.0 && type=="gro") {
377 0 : error("gro files should be in nm");
378 : }
379 17 : if(myunit.getLength()!=1.0 && type=="xtc") {
380 0 : error("xtc files should be in nm");
381 : }
382 17 : if(myunit.getLength()!=1.0 && type=="trr") {
383 0 : error("trr files should be in nm");
384 : }
385 17 : lenunit=getUnits().getLength()/myunit.getLength();
386 227 : } else if(type=="gro" || type=="xtc" || type=="trr") {
387 56 : lenunit=getUnits().getLength();
388 : } else {
389 154 : lenunit=1.0;
390 : }
391 :
392 227 : of.link(*this);
393 227 : of.open(file);
394 : std::string path=of.getPath();
395 227 : log<<" Writing on file "<<path<<"\n";
396 : std::string mode=of.getMode();
397 227 : if(type=="xtc" || type=="trr") {
398 10 : of.close();
399 : // only the first process of each replica writes: the positions are the same
400 : // on all of them, and several processes writing on the same file corrupt it
401 10 : int ok=1;
402 10 : if(comm.Get_rank()==0) {
403 10 : xd=xdrfile::xdrfile_open(path.c_str(),mode.c_str());
404 10 : ok=(xd!=nullptr);
405 : }
406 10 : comm.Bcast(ok,0);
407 10 : if(!ok) {
408 0 : error("cannot open file " + path);
409 : }
410 : }
411 227 : log.printf(" printing the following atoms in %s :", unitname.c_str() );
412 31859 : for(unsigned i=0; i<atoms.size(); ++i) {
413 31632 : log.printf(" %d",atoms[i].serial() );
414 : }
415 227 : log.printf("\n");
416 :
417 227 : if( getNumberOfArguments()>0 ) {
418 40 : if( type!="xyz" ) {
419 0 : error("can only print atomic properties when outputting xyz files");
420 : }
421 :
422 40 : std::vector<std::string> argnames;
423 121 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
424 81 : if( getPntrToArgument(i)->getRank()!=1 ) {
425 0 : error("arguments for xyz output should be vectors");
426 : }
427 81 : if( getPntrToArgument(i)->getNumberOfValues()!=atoms.size() ) {
428 0 : error("number of elements in vector " + getPntrToArgument(i)->getName() + " is not equal to number of atoms output");
429 : }
430 81 : argnames.push_back( getPntrToArgument(i)->getName() );
431 : }
432 40 : std::vector<std::string> str_upper, str_lower;
433 : std::string errors;
434 40 : parseVector("LESS_THAN_OR_EQUAL",str_upper);
435 40 : parseVector("GREATER_THAN_OR_EQUAL",str_lower);
436 40 : if( !bounds.setBounds( getNumberOfArguments(), str_lower, str_upper, errors ) ) {
437 0 : error( errors );
438 : }
439 40 : if( bounds.wereSet() ) {
440 26 : log.printf(" %s \n", bounds.report( argnames ).c_str() );
441 : }
442 40 : }
443 :
444 227 : requestAtoms(atoms, false);
445 227 : auto* infomoldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
446 227 : if( infomoldat ) {
447 59 : log<<" MOLINFO DATA found with label " <<infomoldat->getLabel()<<", using proper atom names\n";
448 59 : names.resize(atoms.size());
449 6563 : for(unsigned i=0; i<atoms.size(); i++)
450 6504 : if(atoms[i].index()<infomoldat->getPDBsize()) {
451 13004 : names[i]=infomoldat->getAtomName(atoms[i]);
452 : }
453 59 : residueNumbers.resize(atoms.size());
454 6563 : for(unsigned i=0; i<residueNumbers.size(); ++i)
455 6504 : if(atoms[i].index()<infomoldat->getPDBsize()) {
456 6502 : residueNumbers[i]=infomoldat->getResidueNumber(atoms[i]);
457 : }
458 59 : residueNames.resize(atoms.size());
459 6563 : for(unsigned i=0; i<residueNames.size(); ++i)
460 6504 : if(atoms[i].index()<infomoldat->getPDBsize()) {
461 13004 : residueNames[i]=infomoldat->getResidueName(atoms[i]);
462 : }
463 : }
464 227 : }
465 :
466 0 : void DumpAtoms::calculateNumericalDerivatives( ActionWithValue* a ) {
467 0 : plumed_merror("this should never be called");
468 : }
469 :
470 9113 : void DumpAtoms::lockRequests() {
471 : ActionWithArguments::lockRequests();
472 : ActionAtomistic::lockRequests();
473 9113 : }
474 :
475 9113 : void DumpAtoms::unlockRequests() {
476 : ActionWithArguments::unlockRequests();
477 : ActionAtomistic::unlockRequests();
478 9113 : }
479 :
480 9018 : void DumpAtoms::update() {
481 9018 : if(type=="xyz") {
482 : unsigned nat=0;
483 8432 : std::vector<double> args( getNumberOfArguments() );
484 77823 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
485 100293 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
486 30902 : args[j] = getPntrToArgument(j)->get(i);
487 : }
488 69391 : if( bounds.check( args ) ) {
489 65585 : nat++;
490 : }
491 : }
492 8432 : of.printf("%d\n",nat);
493 8432 : const Tensor & t(getPbc().getBox());
494 8432 : if(getPbc().isOrthorombic()) {
495 1570 : of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2));
496 : } else {
497 81180 : of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),
498 8118 : lenunit*t(0,0),lenunit*t(0,1),lenunit*t(0,2),
499 8118 : lenunit*t(1,0),lenunit*t(1,1),lenunit*t(1,2),
500 8118 : lenunit*t(2,0),lenunit*t(2,1),lenunit*t(2,2)
501 : );
502 : }
503 25296 : const std::string posFormatString="%s "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz;
504 8432 : const std::string extraFormatString=" "+fmt_xyz;
505 77823 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
506 100293 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
507 30902 : args[j] = getPntrToArgument(j)->get(i);
508 : }
509 69391 : if( !bounds.check(args) ) {
510 3806 : continue;
511 : }
512 : const char* defname="X";
513 : const char* atomName=defname;
514 65585 : if(names.size()>0) {
515 13585 : if(names[i].length()>0) {
516 : atomName=names[i].c_str();
517 : }
518 : }
519 65585 : of.printf(posFormatString.c_str(),
520 : atomName,
521 65585 : lenunit*getPosition(i)(0),
522 65585 : lenunit*getPosition(i)(1),
523 65585 : lenunit*getPosition(i)(2));
524 92681 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
525 27096 : of.printf(extraFormatString.c_str(), getPntrToArgument(j)->get(i) );
526 : }
527 65585 : of.printf("\n");
528 : }
529 9018 : } else if(type=="gro") {
530 466 : const std::string posFormatString="%5d%-5s%5s%5d"+fmt_gro_pos+fmt_gro_pos+fmt_gro_pos+"\n";
531 466 : std::string extraFormatString=" "+fmt_xyz;
532 466 : const Tensor & t(getPbc().getBox());
533 466 : of.printf("Made with PLUMED t=%f\n",getTime()/getUnits().getTime());
534 466 : of.printf("%d\n",getNumberOfAtoms());
535 38404 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
536 : const char* defname="X";
537 : const char* atomName=defname;
538 : int residueNumber=0;
539 37938 : if(names.size()>0)
540 37137 : if(names[i].length()>0) {
541 : atomName=names[i].c_str();
542 : }
543 37938 : if(residueNumbers.size()>0) {
544 37137 : residueNumber=residueNumbers[i];
545 : }
546 37938 : std::string resname="";
547 37938 : if(residueNames.size()>0) {
548 37137 : resname=residueNames[i];
549 : }
550 37938 : of.printf(posFormatString.c_str(),
551 : residueNumber%100000,resname.c_str(),atomName,getAbsoluteIndex(i).serial()%100000,
552 37938 : lenunit*getPosition(i)(0),lenunit*getPosition(i)(1),lenunit*getPosition(i)(2));
553 : }
554 4660 : 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(),
555 466 : lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2),
556 466 : lenunit*t(0,1),lenunit*t(0,2),lenunit*t(1,0),
557 466 : lenunit*t(1,2),lenunit*t(2,0),lenunit*t(2,1));
558 120 : } else if(type=="xtc" || type=="trr") {
559 120 : if(!xd) {
560 0 : return;
561 : }
562 : xdrfile::matrix box;
563 120 : const Tensor & t(getPbc().getBox());
564 120 : int natoms=getNumberOfAtoms();
565 120 : int step=getStep();
566 120 : float time=getTime()/getUnits().getTime();
567 120 : float precision=Tools::fastpow(10.0,iprecision);
568 480 : for(int i=0; i<3; i++)
569 1440 : for(int j=0; j<3; j++) {
570 1080 : box[i][j]=lenunit*t(i,j);
571 : }
572 : // here we cannot use a std::vector<rvec> since it does not compile.
573 : // we thus use a std::unique_ptr<rvec[]>
574 : auto pos = Tools::make_unique<xdrfile::rvec[]>(natoms);
575 6456 : for(int i=0; i<natoms; i++)
576 25344 : for(int j=0; j<3; j++) {
577 19008 : pos[i][j]=lenunit*getPosition(i)(j);
578 : }
579 120 : if(type=="xtc") {
580 72 : write_xtc(xd,natoms,step,time,box,&pos[0],precision);
581 48 : } else if(type=="trr") {
582 48 : write_trr(xd,natoms,step,time,0.0,box,&pos[0],NULL,NULL);
583 : }
584 : } else {
585 0 : plumed_merror("unknown file type "+type);
586 : }
587 : }
588 :
589 454 : DumpAtoms::~DumpAtoms() {
590 227 : if(xd) {
591 10 : xdrfile_close(xd);
592 : }
593 681 : }
594 :
595 :
596 : }
597 : }
|
