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