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