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