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