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/ActionPilot.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/Pbc.h"
26 : #include "tools/File.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/Atoms.h"
29 : #include "tools/Units.h"
30 : #include <cstdio>
31 : #include <memory>
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 : The atoms required are output in a xyz or gro formatted file.
47 : If PLUMED has been compiled with xdrfile support, then also xtc and trr files can be written.
48 : To this aim one should install xdrfile library (http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library).
49 : If the xdrfile library is installed properly the PLUMED configure script should be able to
50 : detect it and enable it.
51 : The type of file is automatically detected from the file extension, but can be also
52 : enforced with TYPE.
53 : Importantly, if your
54 : input file contains actions that edit the atoms position (e.g. \ref WHOLEMOLECULES)
55 : and the DUMPATOMS command appears after this instruction, then the edited
56 : atom positions are output.
57 : You can control the buffering of output using the \ref FLUSH keyword on a separate line.
58 :
59 : Units of the printed file can be controlled with the UNITS keyword. By default PLUMED units as
60 : controlled in the \ref UNITS command are used, but one can override it e.g. with UNITS=A.
61 : Notice that gro/xtc/trr files can only contain coordinates in nm.
62 :
63 : \par Examples
64 :
65 : The following input instructs plumed to print out the positions of atoms
66 : 1-10 together with the position of the center of mass of atoms 11-20 every
67 : 10 steps to a file called file.xyz.
68 : \plumedfile
69 : COM ATOMS=11-20 LABEL=c1
70 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1
71 : \endplumedfile
72 : Notice that the coordinates in the xyz file will be expressed in nm, since these
73 : are the defaults units in PLUMED. If you want the xyz file to be expressed in A, you should use the
74 : following input
75 : \plumedfile
76 : COM ATOMS=11-20 LABEL=c1
77 : DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1 UNITS=A
78 : \endplumedfile
79 : As an alternative, you might want to set all the length used by PLUMED to Angstrom using the \ref UNITS
80 : action. However, this latter choice will affect all your input and output.
81 :
82 : The following input is very similar but dumps a .gro (gromacs) file,
83 : which also contains atom and residue names.
84 : \plumedfile
85 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
86 : # this is required to have proper atom names:
87 : MOLINFO STRUCTURE=reference.pdb
88 : # if omitted, atoms will have "X" name...
89 :
90 : COM ATOMS=11-20 LABEL=c1
91 : DUMPATOMS STRIDE=10 FILE=file.gro ATOMS=1-10,c1
92 : # notice that last atom is a virtual one and will not have
93 : # a correct name in the resulting gro file
94 : \endplumedfile
95 :
96 : The `file.gro` will contain coordinates expressed in nm, since this is the convention for gro files.
97 :
98 : In case you have compiled PLUMED with `xdrfile` library, you might even write xtc or trr files as follows
99 : \plumedfile
100 : COM ATOMS=11-20 LABEL=c1
101 : DUMPATOMS STRIDE=10 FILE=file.xtc ATOMS=1-10,c1
102 : \endplumedfile
103 : Notice that xtc files are significantly smaller than gro and xyz files.
104 :
105 : Finally, consider that gro and xtc file store coordinates with limited precision set by the
106 : `PRECISION` keyword. Default value is 3, which means "3 digits after dot" in nm (1/1000 of a nm).
107 : The following will write a larger xtc file with high resolution coordinates:
108 : \plumedfile
109 : COM ATOMS=11-20 LABEL=c1
110 : DUMPATOMS STRIDE=10 FILE=file.xtc ATOMS=1-10,c1 PRECISION=7
111 : \endplumedfile
112 :
113 :
114 :
115 : */
116 : //+ENDPLUMEDOC
117 :
118 : class DumpAtoms:
119 : public ActionAtomistic,
120 : public ActionPilot {
121 : OFile of;
122 : double lenunit;
123 : int iprecision;
124 : std::vector<std::string> names;
125 : std::vector<unsigned> residueNumbers;
126 : std::vector<std::string> residueNames;
127 : std::string type;
128 : std::string fmt_gro_pos;
129 : std::string fmt_gro_box;
130 : std::string fmt_xyz;
131 : xdrfile::XDRFILE* xd;
132 : public:
133 : explicit DumpAtoms(const ActionOptions&);
134 : ~DumpAtoms();
135 : static void registerKeywords( Keywords& keys );
136 8973 : void calculate() override {}
137 8973 : void apply() override {}
138 : void update() override ;
139 : };
140 :
141 14135 : PLUMED_REGISTER_ACTION(DumpAtoms,"DUMPATOMS")
142 :
143 179 : void DumpAtoms::registerKeywords( Keywords& keys ) {
144 179 : Action::registerKeywords( keys );
145 179 : ActionPilot::registerKeywords( keys );
146 179 : ActionAtomistic::registerKeywords( keys );
147 358 : keys.add("compulsory","STRIDE","1","the frequency with which the atoms should be output");
148 358 : keys.add("atoms", "ATOMS", "the atom indices whose positions you would like to print out");
149 358 : keys.add("compulsory", "FILE", "file on which to output coordinates; extension is automatically detected");
150 358 : keys.add("compulsory", "UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
151 358 : keys.add("optional", "PRECISION","The number of digits in trajectory file");
152 358 : keys.add("optional", "TYPE","file type, either xyz, gro, xtc, or trr, can override an automatically detected file extension");
153 179 : keys.use("RESTART");
154 179 : keys.use("UPDATE_FROM");
155 179 : keys.use("UPDATE_UNTIL");
156 179 : }
157 :
158 175 : DumpAtoms::DumpAtoms(const ActionOptions&ao):
159 : Action(ao),
160 : ActionAtomistic(ao),
161 : ActionPilot(ao),
162 175 : iprecision(3) {
163 : std::vector<AtomNumber> atoms;
164 : std::string file;
165 350 : parse("FILE",file);
166 175 : if(file.length()==0) {
167 0 : error("name out output file was not specified");
168 : }
169 175 : type=Tools::extension(file);
170 175 : log<<" file name "<<file<<"\n";
171 353 : if(type=="gro" || type=="xyz" || type=="xtc" || type=="trr") {
172 157 : log<<" file extension indicates a "<<type<<" file\n";
173 : } else {
174 18 : log<<" file extension not detected, assuming xyz\n";
175 : type="xyz";
176 : }
177 : std::string ntype;
178 350 : parse("TYPE",ntype);
179 175 : if(ntype.length()>0) {
180 5 : if(ntype!="xyz" && ntype!="gro" && ntype!="xtc" && ntype!="trr"
181 : ) {
182 0 : error("TYPE cannot be understood");
183 : }
184 2 : log<<" file type enforced to be "<<ntype<<"\n";
185 : type=ntype;
186 : }
187 :
188 : fmt_gro_pos="%8.3f";
189 : fmt_gro_box="%12.7f";
190 : fmt_xyz="%f";
191 :
192 : std::string precision;
193 350 : parse("PRECISION",precision);
194 175 : if(precision.length()>0) {
195 92 : Tools::convert(precision,iprecision);
196 92 : log<<" with precision "<<iprecision<<"\n";
197 : std::string a,b;
198 92 : Tools::convert(iprecision+5,a);
199 92 : Tools::convert(iprecision,b);
200 184 : fmt_gro_pos="%"+a+"."+b+"f";
201 : fmt_gro_box=fmt_gro_pos;
202 : fmt_xyz=fmt_gro_box;
203 : }
204 :
205 350 : parseAtomList("ATOMS",atoms);
206 :
207 : std::string unitname;
208 350 : parse("UNITS",unitname);
209 175 : if(unitname!="PLUMED") {
210 4 : Units myunit;
211 4 : myunit.setLength(unitname);
212 6 : if(myunit.getLength()!=1.0 && type=="gro") {
213 0 : error("gro files should be in nm");
214 : }
215 6 : if(myunit.getLength()!=1.0 && type=="xtc") {
216 0 : error("xtc files should be in nm");
217 : }
218 6 : if(myunit.getLength()!=1.0 && type=="trr") {
219 0 : error("trr files should be in nm");
220 : }
221 4 : lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
222 419 : } else if(type=="gro" || type=="xtc" || type=="trr") {
223 56 : lenunit=plumed.getAtoms().getUnits().getLength();
224 : } else {
225 115 : lenunit=1.0;
226 : }
227 :
228 175 : checkRead();
229 175 : of.link(*this);
230 175 : of.open(file);
231 : std::string path=of.getPath();
232 175 : log<<" Writing on file "<<path<<"\n";
233 : std::string mode=of.getMode();
234 175 : if(type=="xtc") {
235 6 : of.close();
236 6 : xd=xdrfile::xdrfile_open(path.c_str(),mode.c_str());
237 169 : } else if(type=="trr") {
238 4 : of.close();
239 4 : xd=xdrfile::xdrfile_open(path.c_str(),mode.c_str());
240 : }
241 175 : log.printf(" printing the following atoms in %s :", unitname.c_str() );
242 18836 : for(unsigned i=0; i<atoms.size(); ++i) {
243 18661 : log.printf(" %d",atoms[i].serial() );
244 : }
245 175 : log.printf("\n");
246 175 : requestAtoms(atoms);
247 175 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
248 175 : if( moldat ) {
249 56 : log<<" MOLINFO DATA found with label " <<moldat->getLabel()<<", using proper atom names\n";
250 56 : names.resize(atoms.size());
251 5835 : for(unsigned i=0; i<atoms.size(); i++)
252 5779 : if(atoms[i].index()<moldat->getPDBsize()) {
253 11554 : names[i]=moldat->getAtomName(atoms[i]);
254 : }
255 56 : residueNumbers.resize(atoms.size());
256 5835 : for(unsigned i=0; i<residueNumbers.size(); ++i)
257 5779 : if(atoms[i].index()<moldat->getPDBsize()) {
258 5777 : residueNumbers[i]=moldat->getResidueNumber(atoms[i]);
259 : }
260 56 : residueNames.resize(atoms.size());
261 5835 : for(unsigned i=0; i<residueNames.size(); ++i)
262 5779 : if(atoms[i].index()<moldat->getPDBsize()) {
263 11554 : residueNames[i]=moldat->getResidueName(atoms[i]);
264 : }
265 : }
266 175 : }
267 :
268 8878 : void DumpAtoms::update() {
269 8878 : if(type=="xyz") {
270 8292 : of.printf("%d\n",getNumberOfAtoms());
271 8292 : const Tensor & t(getPbc().getBox());
272 8292 : if(getPbc().isOrthorombic()) {
273 348 : of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2));
274 : } else {
275 16236 : of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),
276 8118 : lenunit*t(0,0),lenunit*t(0,1),lenunit*t(0,2),
277 8118 : lenunit*t(1,0),lenunit*t(1,1),lenunit*t(1,2),
278 8118 : lenunit*t(2,0),lenunit*t(2,1),lenunit*t(2,2)
279 : );
280 : }
281 52086 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
282 : const char* defname="X";
283 : const char* name=defname;
284 43794 : if(names.size()>0)
285 10802 : if(names[i].length()>0) {
286 : name=names[i].c_str();
287 : }
288 87588 : of.printf(("%s "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),name,lenunit*getPosition(i)(0),lenunit*getPosition(i)(1),lenunit*getPosition(i)(2));
289 : }
290 586 : } else if(type=="gro") {
291 466 : const Tensor & t(getPbc().getBox());
292 466 : of.printf("Made with PLUMED t=%f\n",getTime()/plumed.getAtoms().getUnits().getTime());
293 466 : of.printf("%d\n",getNumberOfAtoms());
294 38404 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
295 : const char* defname="X";
296 : const char* name=defname;
297 : unsigned residueNumber=0;
298 37938 : if(names.size()>0)
299 37137 : if(names[i].length()>0) {
300 : name=names[i].c_str();
301 : }
302 37938 : if(residueNumbers.size()>0) {
303 37137 : residueNumber=residueNumbers[i];
304 : }
305 37938 : std::string resname="";
306 37938 : if(residueNames.size()>0) {
307 37137 : resname=residueNames[i];
308 : }
309 75876 : of.printf(("%5u%-5s%5s%5d"+fmt_gro_pos+fmt_gro_pos+fmt_gro_pos+"\n").c_str(),
310 : residueNumber%100000,resname.c_str(),name,getAbsoluteIndex(i).serial()%100000,
311 37938 : lenunit*getPosition(i)(0),lenunit*getPosition(i)(1),lenunit*getPosition(i)(2));
312 : }
313 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(),
314 466 : lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2),
315 466 : lenunit*t(0,1),lenunit*t(0,2),lenunit*t(1,0),
316 466 : lenunit*t(1,2),lenunit*t(2,0),lenunit*t(2,1));
317 168 : } else if(type=="xtc" || type=="trr") {
318 : xdrfile::matrix box;
319 120 : const Tensor & t(getPbc().getBox());
320 120 : int natoms=getNumberOfAtoms();
321 120 : int step=getStep();
322 120 : float time=getTime()/plumed.getAtoms().getUnits().getTime();
323 120 : float precision=Tools::fastpow(10.0,iprecision);
324 480 : for(int i=0; i<3; i++)
325 1440 : for(int j=0; j<3; j++) {
326 1080 : box[i][j]=lenunit*t(i,j);
327 : }
328 : // here we cannot use a std::vector<rvec> since it does not compile.
329 : // we thus use a std::unique_ptr<rvec[]>
330 120 : auto pos = Tools::make_unique<xdrfile::rvec[]>(natoms);
331 6456 : for(int i=0; i<natoms; i++)
332 25344 : for(int j=0; j<3; j++) {
333 19008 : pos[i][j]=lenunit*getPosition(i)(j);
334 : }
335 120 : if(type=="xtc") {
336 72 : write_xtc(xd,natoms,step,time,box,&pos[0],precision);
337 48 : } else if(type=="trr") {
338 48 : write_trr(xd,natoms,step,time,0.0,box,&pos[0],NULL,NULL);
339 : }
340 : } else {
341 0 : plumed_merror("unknown file type "+type);
342 : }
343 8878 : }
344 :
345 350 : DumpAtoms::~DumpAtoms() {
346 175 : if(type=="xtc") {
347 6 : xdrfile_close(xd);
348 169 : } else if(type=="trr") {
349 4 : xdrfile_close(xd);
350 : }
351 700 : }
352 :
353 :
354 : }
355 : }
|
