Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2025 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 "TrajectoryParser.h"
23 :
24 : #include "FileBase.h"
25 : #include "IFile.h"
26 : #include "Tools.h"
27 : #include "xdrfile/xdrfile.h"
28 : #include "xdrfile/xdrfile_trr.h"
29 : #include "xdrfile/xdrfile_xtc.h"
30 :
31 : // when using molfile plugin
32 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
33 : #ifndef __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS
34 : /* Use the internal ones. Alternatively:
35 : * ifeq (,$(findstring __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS,$(CPPFLAGS)))
36 : * CPPFLAGS+=-I../molfile
37 : */
38 : #include "molfile/libmolfile_plugin.h"
39 : #include "molfile/molfile_plugin.h"
40 : using namespace PLMD::molfile;
41 : #else
42 : #include <libmolfile_plugin.h>
43 : #include <molfile_plugin.h>
44 : #endif
45 :
46 : static int register_cb(void *v, vmdplugin_t *p);
47 : struct molFilePlugins {
48 : std::vector<molfile_plugin_t *> plugins{};
49 : std::map <std::string, unsigned> pluginmap{};
50 1084768 : static molFilePlugins& get() {
51 1084768 : static molFilePlugins me;
52 1084768 : return me;
53 : }
54 1276 : std::vector<std::string> getMolfilePluginsnames() {
55 1276 : std::vector<std::string> toRet(plugins.size());
56 12760 : for(unsigned i=0; i<plugins.size(); i++) {
57 22968 : toRet[i] = std::string(plugins[i]->name);
58 : }
59 1276 : return toRet;
60 0 : }
61 : private:
62 : molFilePlugins() {};
63 : };
64 197676 : static int register_cb(void *v, vmdplugin_t *p) {
65 395352 : const auto ret = molFilePlugins::get().pluginmap.insert ( std::pair<std::string,unsigned>(
66 197676 : std::string(p->name),
67 197676 : molFilePlugins::get().plugins.size()) );
68 :
69 197676 : if (ret.second) {
70 65889 : molFilePlugins::get().plugins.push_back(reinterpret_cast<molfile_plugin_t *>(p));
71 : #ifdef PLUMED_MOLFILE_PLUGIN_DEBUG
72 : std::cerr<<"MOLFILE: loading plugin "<<p->name<<" number "<<molFilePlugins::get().plugins.size()-1<<std::endl;
73 : } else {
74 : std::cerr<<"MOLFILE: found duplicate plugin for "<<p->name<<" : not inserted "<<std::endl;
75 : #endif
76 : }
77 197676 : return VMDPLUGIN_SUCCESS;
78 : }
79 : #endif
80 :
81 :
82 : namespace PLMD {
83 :
84 1265 : class fileParser {
85 : protected:
86 : int natoms {-1};
87 : public:
88 : virtual ~fileParser()=default;
89 : int nOfAtoms() const {
90 78791 : return natoms;
91 : }
92 : virtual std::optional<std::string> init(std::string_view fmt,
93 : std::string_view fname,
94 : int command_line_natoms=-1)=0;
95 : virtual std::optional<std::string> init(FILE* /*fileHandle*/) {
96 : return "this trajectory format cannot be read from standard input";
97 : };
98 :
99 : virtual std::optional<std::string> readHeader(
100 : long long int &step,
101 : double &timeStep
102 : )=0;
103 :
104 : virtual std::optional<std::string> readAtoms( int stride,
105 : bool dont_read_pbc,
106 : bool debug_pd,
107 : int pd_start,
108 : int pd_nlocal,
109 : long long int &step,
110 : double* masses,
111 : double* charges,
112 : double* coordinates,
113 : double* cell )=0;
114 :
115 : virtual std::optional<std::string> readHeader(
116 : long long int &step,
117 : float &timeStep
118 : ) =0;
119 : virtual std::optional<std::string> readAtoms(int stride,
120 : bool dont_read_pbc,
121 : bool debug_pd,
122 : int pd_start,
123 : int pd_nlocal,
124 : long long int &step,
125 : float* masses,
126 : float* charges,
127 : float* coordinates,
128 : float* cell ) =0;
129 :
130 0 : virtual std::optional<std::string> rewind() {
131 0 : return "rewind not (yet) implemented for this kind of trajectory file";
132 : }
133 :
134 : };
135 : namespace {
136 :
137 : #define READATOMS \
138 : std::optional<std::string> readAtoms(int stride,bool dont_read_pbc, \
139 : bool debug_pd,int pd_start,int pd_nlocal,long long int &step, \
140 : double* masses,double* charges,double* coordinates,double* cell ) override { \
141 : return readAtoms_t<double>(stride,dont_read_pbc,debug_pd,pd_start,pd_nlocal, \
142 : step,masses,charges,coordinates,cell ); \
143 : } \
144 : \
145 : std::optional<std::string> readAtoms(int stride,bool dont_read_pbc, \
146 : bool debug_pd,int pd_start,int pd_nlocal,long long int &step, \
147 : float* masses,float* charges,float* coordinates,float* cell ) override { \
148 : return readAtoms_t<float>(stride,dont_read_pbc,debug_pd,pd_start,pd_nlocal, \
149 : step,masses,charges,coordinates,cell ); \
150 : }
151 : struct file_deleter {
152 : void operator()(FILE* fp) {
153 0 : std::fclose(fp);
154 0 : }
155 : };
156 :
157 : using safeFile=std::unique_ptr<std::FILE, file_deleter>;
158 748 : class xyzParser final :public fileParser {
159 : // not owning pointer for "-"mode
160 : std::FILE* fp{nullptr};
161 : // owning pointer otherwise
162 : safeFile fp_owned{nullptr,{}};
163 :
164 : template <typename real>
165 41359 : std::optional<std::string> readAtoms_t(int stride,
166 : bool dont_read_pbc,
167 : bool debug_pd,
168 : int pd_start,
169 : int pd_nlocal,
170 : long long int &step,
171 : real* masses,
172 : real* charges,
173 : real* coordinates,
174 : real* cell ) {
175 : std::string line;
176 : //cell
177 41359 : if(!Tools::getline(fp,line)) {
178 : return "premature end of trajectory file";
179 : }
180 41359 : std::vector<double> celld(9,0.0);
181 41359 : if(!dont_read_pbc) {
182 40975 : std::vector<std::string> words=Tools::getWords(line);
183 40975 : if(words.size()==3) {
184 40277 : Tools::convert(words[0],celld[0]);
185 40277 : Tools::convert(words[1],celld[4]);
186 40277 : Tools::convert(words[2],celld[8]);
187 698 : } else if(words.size()==9) {
188 698 : Tools::convert(words[0],celld[0]);
189 698 : Tools::convert(words[1],celld[1]);
190 698 : Tools::convert(words[2],celld[2]);
191 698 : Tools::convert(words[3],celld[3]);
192 698 : Tools::convert(words[4],celld[4]);
193 698 : Tools::convert(words[5],celld[5]);
194 698 : Tools::convert(words[6],celld[6]);
195 698 : Tools::convert(words[7],celld[7]);
196 698 : Tools::convert(words[8],celld[8]);
197 : } else {
198 : return "needed box in second line of xyz file";
199 : }
200 :
201 40975 : }
202 : std::copy(celld.begin(),celld.end(),cell);
203 : //atoms, finally
204 2353638 : for(int i=0; i<natoms; i++) {
205 : double cc[3];
206 2312279 : if(!Tools::getline(fp,line)) {
207 0 : return "premature end of trajectory file";
208 : }
209 : char dummy[1000];
210 2312279 : int ret=std::sscanf(line.c_str(),"%999s %100lf %100lf %100lf",
211 : dummy,&cc[0],&cc[1],&cc[2]);
212 2312279 : if(ret!=4) {
213 0 : return "cannot read line"+line;
214 : }
215 2312279 : if(!debug_pd || ( i>=pd_start && i<pd_start+pd_nlocal) ) {
216 2292839 : coordinates[3*i]=real(cc[0]);
217 2292839 : coordinates[3*i+1]=real(cc[1]);
218 2292839 : coordinates[3*i+2]=real(cc[2]);
219 : }
220 : }
221 41359 : return std::nullopt;
222 41359 : }
223 : public:
224 748 : std::optional<std::string> init(std::string_view fmt,
225 : std::string_view fname,
226 : int command_line_natoms=-1) override {
227 748 : std::string trajectoryFile=std::string(fname);
228 748 : fp_owned.reset(std::fopen(trajectoryFile.c_str(),"r"));
229 748 : fp = fp_owned.get();
230 748 : if(!fp) {
231 2 : return "ERROR: Error opening trajectory file "+trajectoryFile;
232 : }
233 747 : return std::nullopt;
234 : }
235 0 : std::optional<std::string> init(FILE* fileHandle) override {
236 0 : fp = fileHandle;
237 0 : return std::nullopt;
238 : }
239 42100 : std::optional<std::string> readHeader(
240 : long long int &/*step*/,
241 : double &/*timeStep*/
242 : ) override {
243 : std::string line;
244 42100 : if(!Tools::getline(fp,line)) {
245 : return "EOF";
246 : }
247 41409 : std::sscanf(line.c_str(),"%100d",&natoms);
248 41409 : return std::nullopt;
249 : }
250 :
251 5 : std::optional<std::string> readHeader(
252 : long long int &step,
253 : float &/*timeStep*/
254 : ) override {
255 : double ts;
256 : //ignore this -Wmaybe-uninitialized
257 5 : auto msg=readHeader(step,ts);
258 5 : return msg;
259 : }
260 : //see the template readAtoms_t for the implementation
261 41359 : READATOMS;
262 1 : std::optional<std::string> rewind() override {
263 1 : int error = std::fseek(fp,0,SEEK_SET);
264 1 : if (error) {
265 : return "ERROR: Error rewinding trajectory file";
266 : }
267 1 : return std::nullopt;
268 : }
269 : };
270 :
271 3 : class dlp4Parser final:public fileParser {
272 : // not owning pointer for "-"mode
273 : std::FILE* fp{nullptr};
274 : // owning pointer otherwise
275 : safeFile fp_owned{nullptr,{}};
276 : int lvl{0};
277 : template <typename real>
278 30 : std::optional<std::string> readAtoms_t(int stride,
279 : bool dont_read_pbc,
280 : bool debug_pd,
281 : int pd_start,
282 : int pd_nlocal,
283 : long long int &step,
284 : real* masses,
285 : real* charges,
286 : real* coordinates,
287 : real* cell ) {
288 : std::string line;
289 :
290 : //cell
291 30 : std::vector<double> celld(9,0.0);
292 30 : if(!dont_read_pbc) {
293 30 : if(!Tools::getline(fp,line)) {
294 : return "error reading vector a of cell";
295 : }
296 30 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[0],&celld[1],&celld[2]);
297 30 : if(!Tools::getline(fp,line)) {
298 : return "error reading vector b of cell";
299 : }
300 30 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[3],&celld[4],&celld[5]);
301 30 : if(!Tools::getline(fp,line)) {
302 : return "error reading vector c of cell";
303 : }
304 30 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[6],&celld[7],&celld[8]);
305 : std::transform(celld.begin(),celld.end(),celld.begin(),
306 : [](double x) {
307 : //converting to nm
308 270 : return x/10.0;
309 : });
310 : std::copy(celld.begin(),celld.end(),cell);
311 : }
312 : // reading the frame
313 990 : for(int i=0; i<natoms; i++) {
314 : double cc[3];
315 960 : bool ok=Tools::getline(fp,line);
316 960 : if(!ok) {
317 0 : return "premature end of trajectory file";
318 : }
319 :
320 : char dummy[9];
321 : int idummy;
322 : double m,c;
323 960 : auto res=std::sscanf(line.c_str(),"%8s %d %lf %lf",dummy,&idummy,&m,&c);
324 960 : if (res != 4) {
325 0 : return "cannot read line"+line;
326 : }
327 960 : masses[i]=real(m);
328 960 : charges[i]=real(c);
329 960 : if(!Tools::getline(fp,line)) {
330 : return "error reading coordinates";
331 : }
332 960 : std::sscanf(line.c_str(),"%lf %lf %lf",&cc[0],&cc[1],&cc[2]);
333 960 : cc[0]*=0.1;
334 960 : cc[1]*=0.1;
335 960 : cc[2]*=0.1;
336 960 : if(lvl>0) {
337 960 : if(!Tools::getline(fp,line)) {
338 : return "error skipping velocities";
339 : }
340 : }
341 960 : if(lvl>1) {
342 960 : if(!Tools::getline(fp,line)) {
343 : return "error skipping forces";
344 : }
345 : }
346 960 : if(!debug_pd || ( i>=pd_start && i<pd_start+pd_nlocal) ) {
347 960 : coordinates[3*i]=real(cc[0]);
348 960 : coordinates[3*i+1]=real(cc[1]);
349 960 : coordinates[3*i+2]=real(cc[2]);
350 : }
351 : }
352 :
353 30 : return std::nullopt;
354 30 : }
355 :
356 : public:
357 3 : std::optional<std::string> init(std::string_view fmt,
358 : std::string_view fname,
359 : int command_line_natoms=-1) override {
360 3 : std::string trajectoryFile=std::string(fname);
361 3 : fp_owned.reset(std::fopen(trajectoryFile.c_str(),"r"));
362 3 : fp = fp_owned.get();
363 3 : if(!fp) {
364 0 : return "ERROR: Error opening trajectory file "+trajectoryFile;
365 : }
366 : // dlp4 header
367 : std::string line;
368 3 : if(!Tools::getline(fp,line)) {
369 : return "error reading title";
370 : }
371 3 : if(!Tools::getline(fp,line)) {
372 : return "error reading atoms";
373 : }
374 : // in the original driver this pb is only read here, and then not used anymore
375 3 : int pb=1;
376 3 : std::sscanf(line.c_str(),"%d %d %d",&lvl,&pb,&natoms);
377 3 : return std::nullopt;
378 : }
379 0 : std::optional<std::string> init(FILE* fileHandle) override {
380 0 : fp = fileHandle;
381 0 : return std::nullopt;
382 : }
383 :
384 33 : std::optional<std::string> readHeader (
385 : long long int &step,
386 : double &timeStep
387 : ) override {
388 : std::string line;
389 33 : if(!Tools::getline(fp,line)) {
390 : return "EOF";
391 : }
392 : char xa[9];
393 : int xb,xc,xd;
394 : double t;
395 30 : auto res = std::sscanf(line.c_str(),"%8s %lld %d %d %d %lf",xa,&step,&xb,&xc,&xd,&t);
396 30 : if (res!=6) {
397 : return "error reading header";
398 : }
399 30 : timeStep=t;
400 30 : return std::nullopt;
401 : }
402 :
403 0 : std::optional<std::string> readHeader(
404 : long long int &step,
405 : float &timeStep
406 : ) override {
407 0 : double ts=timeStep;
408 0 : auto msg=readHeader(step,ts);
409 0 : timeStep = float(ts);
410 0 : return msg;
411 : }
412 : //see the template readAtoms_t for the implementation
413 30 : READATOMS;
414 0 : std::optional<std::string> rewind() override {
415 0 : int error = std::fseek(fp,0,SEEK_SET);
416 0 : if (error) {
417 : return "ERROR: Error rewinding trajectory file";
418 : }
419 0 : return std::nullopt;
420 : }
421 : };
422 :
423 137 : class groParser final:public fileParser {
424 : // not owning pointer for "-"mode
425 : std::FILE* fp{nullptr};
426 : // owning pointer otherwise
427 : safeFile fp_owned{nullptr,{}};
428 : template <typename real>
429 12688 : std::optional<std::string> readAtoms_t(int stride,
430 : bool dont_read_pbc,
431 : bool debug_pd,
432 : int pd_start,
433 : int pd_nlocal,
434 : long long int &step,
435 : real* masses,
436 : real* charges,
437 : real* coordinates,
438 : real* cell ) {
439 : std::string line;
440 : int ddist=0;
441 911369 : for(int i=0; i<natoms; i++) {
442 898681 : if(!Tools::getline(fp,line)) {
443 0 : return "premature end of trajectory file";
444 : }
445 : double cc[3];
446 : // do the gromacs way
447 898681 : if(!i) {
448 : //
449 : // calculate the distance between dots (as in gromacs gmxlib/confio.c, routine get_w_conf )
450 : //
451 : const char *p1, *p2, *p3;
452 : p1 = std::strchr(line.c_str(), '.');
453 12688 : if (p1 == NULL) {
454 : return "seems there are no coordinates in the gro file";
455 : }
456 12688 : p2 = std::strchr(&p1[1], '.');
457 12688 : if (p2 == NULL) {
458 : return "seems there is only one coordinates in the gro file";
459 : }
460 12688 : ddist = p2 - p1;
461 12688 : p3 = std::strchr(&p2[1], '.');
462 12688 : if (p3 == NULL) {
463 : return "seems there are only two coordinates in the gro file";
464 : }
465 12688 : if (p3 - p2 != ddist) {
466 : return "not uniform spacing in fields in the gro file";
467 : }
468 : }
469 898681 : Tools::convert(line.substr(20,ddist),cc[0]);
470 898681 : Tools::convert(line.substr(20+ddist,ddist),cc[1]);
471 898681 : Tools::convert(line.substr(20+ddist+ddist,ddist),cc[2]);
472 898681 : if(!debug_pd || ( i>=pd_start && i<pd_start+pd_nlocal) ) {
473 898681 : coordinates[3*i]=real(cc[0]);
474 898681 : coordinates[3*i+1]=real(cc[1]);
475 898681 : coordinates[3*i+2]=real(cc[2]);
476 : }
477 :
478 : }
479 12688 : if(!Tools::getline(fp,line)) {
480 : return "premature end of trajectory file";
481 : }
482 12688 : std::vector<std::string> words=Tools::getWords(line);
483 12688 : if(words.size()<3) {
484 : return "cannot understand box format";
485 : }
486 12688 : Tools::convert(words[0],cell[0]);
487 12688 : Tools::convert(words[1],cell[4]);
488 12688 : Tools::convert(words[2],cell[8]);
489 12688 : if(words.size()>3) {
490 531 : Tools::convert(words[3],cell[1]);
491 : }
492 12688 : if(words.size()>4) {
493 531 : Tools::convert(words[4],cell[2]);
494 : }
495 12688 : if(words.size()>5) {
496 531 : Tools::convert(words[5],cell[3]);
497 : }
498 12688 : if(words.size()>6) {
499 531 : Tools::convert(words[6],cell[5]);
500 : }
501 12688 : if(words.size()>7) {
502 531 : Tools::convert(words[7],cell[6]);
503 : }
504 12688 : if(words.size()>8) {
505 531 : Tools::convert(words[8],cell[7]);
506 : }
507 12688 : return std::nullopt;
508 12688 : }
509 :
510 : public:
511 137 : std::optional<std::string> init(std::string_view fmt,
512 : std::string_view fname,
513 : int command_line_natoms=-1) override {
514 137 : std::string trajectoryFile=std::string(fname);
515 137 : fp_owned.reset(std::fopen(trajectoryFile.c_str(),"r"));
516 137 : fp = fp_owned.get();
517 137 : if(!fp) {
518 0 : return "ERROR: Error opening trajectory file "+trajectoryFile;
519 : }
520 137 : return std::nullopt;
521 : }
522 0 : std::optional<std::string> init(FILE* fileHandle) override {
523 0 : fp = fileHandle;
524 0 : return std::nullopt;
525 : }
526 :
527 12822 : std::optional<std::string> readHeader (
528 : long long int &/**/,
529 : double &/**/
530 : ) override {
531 : std::string line;
532 12822 : if(!Tools::getline(fp,line)) {
533 : return "EOF";
534 : }
535 12688 : if(!Tools::getline(fp,line)) {
536 : return "premature end of trajectory file";
537 : }
538 12688 : std::sscanf(line.c_str(),"%100d",&natoms);
539 :
540 12688 : return std::nullopt;
541 : }
542 :
543 0 : std::optional<std::string> readHeader(
544 : long long int &step,
545 : float &/**/
546 : ) override {
547 : double ts;
548 0 : auto msg=readHeader(step,ts);
549 0 : return msg;
550 : }
551 : //see the template readAtoms_t for the implementation
552 12688 : READATOMS;
553 0 : std::optional<std::string> rewind() override {
554 0 : int error = std::fseek(fp,0,SEEK_SET);
555 0 : if (error) {
556 : return "ERROR: Error rewinding trajectory file";
557 : }
558 0 : return std::nullopt;
559 : }
560 : };
561 :
562 :
563 : struct xd_deleter {
564 : void operator()(PLMD::xdrfile::XDRFILE* xd) {
565 0 : xdrfile::xdrfile_close(xd);
566 0 : }
567 : };
568 :
569 : enum class xdType {trr,xtc};
570 : template <xdType is>
571 14 : class xdParser final: public fileParser {
572 : // owning pointer
573 : std::unique_ptr<PLMD::xdrfile::XDRFILE, xd_deleter> xd{nullptr,{}};
574 :
575 : template <typename real>
576 133 : std::optional<std::string> readAtoms_t(int stride,
577 : bool dont_read_pbc,
578 : bool debug_pd,
579 : int pd_start,
580 : int pd_nlocal,
581 : long long int &step,
582 : real* masses,
583 : real* charges,
584 : real* coordinates,
585 : real* cell ) {
586 : int localstep;
587 : float time;
588 : xdrfile::matrix box;
589 : // here we cannot use a std::vector<rvec> since it does not compile.
590 : // we thus use a std::unique_ptr<rvec[]>
591 : auto pos=Tools::make_unique<xdrfile::rvec[]>(natoms);
592 : int ret=xdrfile::exdrOK;
593 : if constexpr(is == xdType::xtc) {
594 : float prec;
595 52 : ret=xdrfile::read_xtc(xd.get(),natoms,&localstep,&time,box,pos.get(),&prec);
596 : }
597 : if constexpr(is == xdType::trr) {
598 : float lambda;
599 81 : ret=xdrfile::read_trr(xd.get(),natoms,&localstep,&time,&lambda,box,pos.get(),NULL,NULL);
600 : }
601 133 : if(stride==0) {
602 30 : step=localstep;
603 : }
604 133 : if(ret==xdrfile::exdrENDOFFILE) {
605 : return "EOF";
606 : }
607 130 : if(ret!=xdrfile::exdrOK) {
608 : return "EOF";
609 : }
610 480 : for(unsigned i=0; i<3; i++) {
611 1440 : for(unsigned j=0; j<3; j++) {
612 1080 : cell[3*i+j]=box[i][j];
613 : }
614 : }
615 : //xdrfiles are not used with dd ?
616 7825 : for(int i=0; i<natoms; i++) {
617 30820 : for(unsigned j=0; j<3; j++) {
618 23115 : coordinates[3*i+j]=real(pos[i][j]);
619 : }
620 : }
621 120 : return std::nullopt;
622 : }
623 :
624 : public:
625 13 : std::optional<std::string> init(std::string_view fmt,
626 : std::string_view fname,
627 : int command_line_natoms=-1) override {
628 : using trajfmt=TrajectoryParser::trajfmt;
629 13 : std::string trajectoryFile=std::string(fname);
630 13 : auto trajectory_fmt =TrajectoryParser::FMTfromString(fmt);
631 13 : xd.reset(xdrfile::xdrfile_open(trajectoryFile.c_str(),"r"));
632 13 : if(!xd) {
633 0 : return "ERROR: Error opening trajectory file "+trajectoryFile;
634 : }
635 : if constexpr(is == xdType::xtc) {
636 3 : plumed_assert(trajectory_fmt==trajfmt::xdr_xtc) << "trajectory type should be xtc";
637 3 : xdrfile::read_xtc_natoms(&trajectoryFile[0],&natoms);
638 : }
639 : if constexpr(is == xdType::trr) {
640 10 : plumed_assert(trajectory_fmt==trajfmt::xdr_trr) << "trajectory type should be trr";
641 10 : xdrfile::read_trr_natoms(&trajectoryFile[0],&natoms);
642 : }
643 13 : return std::nullopt;
644 : }
645 :
646 1 : std::optional<std::string> init(FILE* /*fileHandle*/) override {
647 2 : plumed_assert(false) << "xd readers are not compatible with reading from stdin";
648 : return std::nullopt;
649 : }
650 :
651 133 : std::optional<std::string> readHeader (long long int &/*step*/,double &/*timeStep*/) override {
652 133 : return std::nullopt;
653 : }
654 :
655 0 : std::optional<std::string> readHeader(long long int &/*step*/,float &/*timeStep*/) override {
656 0 : return std::nullopt;
657 : }
658 : //see the template readAtoms_t for the implementation
659 133 : READATOMS;
660 : };
661 :
662 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
663 : class molfileParser final: public fileParser {
664 :
665 : molfile::molfile_plugin_t *api=NULL;
666 :
667 : molfile_timestep_t ts_in;
668 : std::string trajectoryFile;
669 : std::string trajectory_fmt_str;
670 : //Initializing ts_in_coords should prevent that ts_in_coords.data() results in UB (and passes codecheck)
671 : std::vector<float> ts_in_coords{0.0f,0.0f,0.0f};
672 : struct molfile_deleter {
673 : molfile::molfile_plugin_t *api=NULL;
674 362 : void operator()(void* h) {
675 362 : if(h) {
676 362 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
677 362 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) {
678 498 : lck=Tools::molfile_lock();
679 : }
680 362 : api->close_file_read(h);
681 362 : }
682 362 : }
683 : };
684 : std::unique_ptr<void, molfile_deleter> h_in{nullptr,{nullptr}};
685 :
686 :
687 : template <typename real>
688 23410 : std::optional<std::string> readAtoms_t(int stride,
689 : bool dont_read_pbc,
690 : bool debug_pd,
691 : int pd_start,
692 : int pd_nlocal,
693 : long long int &step,
694 : real* masses,
695 : real* charges,
696 : real* coordinates,
697 : real* cell ) {
698 23410 : if(!dont_read_pbc) {
699 : //ts_in has been read in the previous step
700 23367 : if(ts_in.A>0.0) { // this is negative if molfile does not provide box
701 : // info on the cell: convert using pbcset.tcl from pbctools in vmd distribution
702 : constexpr real r180=real(180.);
703 23329 : real cosBC=cos(real(ts_in.alpha)*pi/r180);
704 23329 : real cosAC=std::cos(real(ts_in.beta)*pi/r180);
705 23329 : real cosAB=std::cos(real(ts_in.gamma)*pi/r180);
706 23329 : real sinAB=std::sin(real(ts_in.gamma)*pi/r180);
707 23329 : real Ax=real(ts_in.A);
708 23329 : real Bx=real(ts_in.B)*cosAB;
709 23329 : real By=real(ts_in.B)*sinAB;
710 23329 : real Cx=real(ts_in.C)*cosAC;
711 23329 : real Cy=(real(ts_in.C)*real(ts_in.B)*cosBC-Cx*Bx)/By;
712 23329 : real Cz=std::sqrt(real(ts_in.C)*real(ts_in.C)-Cx*Cx-Cy*Cy);
713 23329 : cell[0]=Ax/10.;
714 23329 : cell[1]=0.;
715 23329 : cell[2]=0.;
716 23329 : cell[3]=Bx/10.;
717 23329 : cell[4]=By/10.;
718 23329 : cell[5]=0.;
719 23329 : cell[6]=Cx/10.;
720 23329 : cell[7]=Cy/10.;
721 23329 : cell[8]=Cz/10.;
722 : } else {
723 38 : cell[0]=0.0;
724 38 : cell[1]=0.0;
725 38 : cell[2]=0.0;
726 38 : cell[3]=0.0;
727 38 : cell[4]=0.0;
728 38 : cell[5]=0.0;
729 38 : cell[6]=0.0;
730 38 : cell[7]=0.0;
731 38 : cell[8]=0.0;
732 : }
733 : }
734 : // info on coords
735 : // the order is xyzxyz...
736 18638785 : for(int i=0; i<3*natoms; i++) {
737 18615375 : coordinates[i]=real(ts_in.coords[i])/real(10.); //convert to nm
738 : //cerr<<"COOR "<<coordinates[i]<<endl;
739 : }
740 23410 : return std::nullopt;
741 : }
742 :
743 : public:
744 363 : molfileParser() {
745 363 : ts_in.coords=NULL;
746 363 : ts_in.velocities=NULL;
747 363 : ts_in.A=-1; // we use this to check whether cell is provided or not
748 363 : }
749 363 : std::optional<std::string> init(std::string_view fmt,
750 : std::string_view fname,
751 : int command_line_natoms=-1) override {
752 726 : trajectoryFile = std::string(fname);
753 363 : ts_in.coords=ts_in_coords.data();
754 363 : ts_in.velocities=NULL;
755 363 : ts_in.A=-1; // we use this to check whether cell is provided or not
756 363 : trajectory_fmt_str=std::string(fmt);
757 : bool found=false;
758 3630 : for(unsigned i=0; i<molFilePlugins::get().plugins.size(); ++i) {
759 3267 : if(fmt == molFilePlugins::get().plugins[i]->name) {
760 : found=true;
761 362 : api = molFilePlugins::get().plugins[i];
762 : }
763 : }
764 363 : if(!found) {
765 2 : return "trajectory format "+trajectory_fmt_str+" is not supported by the molfile plugin";
766 : }
767 362 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
768 362 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) {
769 498 : lck=Tools::molfile_lock();
770 : }
771 362 : h_in = std::unique_ptr<void, molfile_deleter> {
772 362 : api->open_file_read(trajectoryFile.c_str(), trajectory_fmt_str.c_str(), &natoms),
773 362 : {api}
774 : };
775 362 : if(natoms==MOLFILE_NUMATOMS_UNKNOWN) {
776 3 : if(command_line_natoms>=0) {
777 3 : natoms=command_line_natoms;
778 : } else {
779 : return "this file format does not provide number of atoms; use --natoms on the command line";
780 : }
781 : }
782 362 : ts_in_coords.resize(3*natoms);
783 362 : ts_in.coords = ts_in_coords.data();
784 :
785 362 : return std::nullopt;
786 362 : }
787 :
788 0 : std::optional<std::string> init(FILE* /*fileHandle*/) override {
789 0 : plumed_assert(false) << "molfile readers are not compatible with reading from stdin";
790 : return std::nullopt;
791 : }
792 :
793 23772 : std::optional<std::string> readHeader (long long int &/*step*/,double &/*timeStep*/) override {
794 23772 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
795 23772 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) {
796 46860 : lck=Tools::molfile_lock();
797 : }
798 : int rc;
799 23772 : rc = api->read_next_timestep(h_in.get(), natoms, &ts_in);
800 23772 : if(rc==MOLFILE_EOF) {
801 : return "EOF";
802 : }
803 23410 : return std::nullopt;
804 23772 : }
805 :
806 0 : std::optional<std::string> readHeader(long long int &/*step*/,float &/*timeStep*/) override {
807 0 : return std::nullopt;
808 : }
809 : //see the template readAtoms_t for the implementation
810 23410 : READATOMS;
811 : };
812 : #endif
813 :
814 : #undef READATOMS
815 :
816 1265 : std::unique_ptr<fileParser> parserFromFMT(TrajectoryParser::trajfmt fmt) {
817 : using trajfmt = TrajectoryParser::trajfmt;
818 1265 : switch (fmt) {
819 363 : case trajfmt::molfile: {
820 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
821 363 : return std::make_unique<molfileParser>(molfileParser {});
822 : #else
823 : plumed_assert(false) <<"Plumed has not been compiled with molfile,\n"
824 : "it is impossible to read a file with the molfile plugin";
825 : #endif
826 : }
827 4 : case trajfmt::xdr_xtc: {
828 4 : return std::make_unique<xdParser<xdType::xtc>>();
829 : }
830 10 : case trajfmt::xdr_trr: {
831 10 : return std::make_unique<xdParser<xdType::trr>>();
832 : }
833 748 : case trajfmt::xyz: {
834 748 : return std::make_unique<xyzParser>();
835 : }
836 137 : case trajfmt::gro: {
837 137 : return std::make_unique<groParser>();
838 : }
839 3 : case trajfmt::dlp4: {
840 3 : return std::make_unique<dlp4Parser>();
841 : }
842 0 : default: {
843 0 : plumed_assert(false) <<"you should not be here" ;
844 : }
845 : }
846 : }
847 : }
848 :
849 21964 : void TrajectoryParser::registerKeywords(Keywords& keys) {
850 21964 : keys.add("atoms","--ixyz","the trajectory in xyz format");
851 21964 : keys.add("atoms","--igro","the trajectory in gro format");
852 21964 : keys.add("atoms","--idlp4","the trajectory in DL_POLY_4 format");
853 21964 : keys.add("atoms","--ixtc","the trajectory in xtc format (xdrfile implementation)");
854 21964 : keys.add("atoms","--itrr","the trajectory in trr format (xdrfile implementation)");
855 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
856 21964 : MOLFILE_INIT_ALL
857 21964 : MOLFILE_REGISTER_ALL(NULL, register_cb)
858 219640 : for(unsigned i=0; i<molFilePlugins::get().plugins.size(); i++) {
859 197676 : std::string kk="--mf_"+std::string(molFilePlugins::get().plugins[i]->name);
860 395352 : std::string mm=" molfile: the trajectory in "+std::string(molFilePlugins::get().plugins[i]->name)+" format " ;
861 197676 : keys.add("atoms",kk,mm);
862 : }
863 : #endif
864 21964 : }
865 :
866 1341 : TrajectoryParser::TrajectoryParser()=default;
867 1341 : TrajectoryParser::~TrajectoryParser()=default;
868 1276 : std::vector<std::string> TrajectoryParser::getMolfilePluginsnames() {
869 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
870 1276 : return molFilePlugins::get().getMolfilePluginsnames();
871 : #else
872 : return {};
873 : #endif
874 : }
875 1276 : std::vector<std::string> TrajectoryParser::trajectoryOptions() {
876 1276 : return {"xyz","gro","dlp4","xtc","trr"};
877 : }
878 :
879 :
880 1278 : TrajectoryParser::trajfmt TrajectoryParser::FMTfromString(std::string_view fmt) {
881 1278 : if(fmt=="xyz") {
882 : return trajfmt::xyz;
883 : }
884 529 : if(fmt=="gro") {
885 : return trajfmt::gro;
886 : }
887 392 : if(fmt=="dlp4") {
888 : return trajfmt::dlp4;
889 : }
890 389 : if(fmt=="xtc") {
891 : return trajfmt::xdr_xtc;
892 : }
893 138 : if(fmt=="trr") {
894 25 : return trajfmt::xdr_trr;
895 : }
896 : return trajfmt::error;
897 : }
898 0 : std::string TrajectoryParser::toString(trajfmt fmt) {
899 0 : switch (fmt) {
900 : case trajfmt::molfile: {
901 0 : return "molfile";
902 : }
903 : case trajfmt::xdr_xtc: {
904 0 : return "xtc";
905 : }
906 : case trajfmt::xdr_trr: {
907 0 : return "trr";
908 : }
909 : case trajfmt::xyz: {
910 0 : return "xyz";
911 : }
912 : case trajfmt::gro: {
913 0 : return "gro";
914 : }
915 : case trajfmt::dlp4: {
916 0 : return "dlp4";
917 : }
918 : case trajfmt::error: {
919 0 : return "error";
920 : }
921 : }
922 :
923 0 : return "unknown";
924 : }
925 :
926 1 : std::optional<std::string> TrajectoryParser::init(std::string_view fmt, FILE* in) {
927 1 : auto trajectory_fmt =FMTfromString(fmt);
928 1 : if (trajectory_fmt==trajfmt::error) {
929 : return "Error: wrong input format specified";
930 : }
931 2 : parser = parserFromFMT(trajectory_fmt);
932 : //not owning, will not be deleted on destruction
933 1 : return parser->init(in);
934 : }
935 1264 : std::optional<std::string> TrajectoryParser::init(std::string_view fmt,
936 : std::string_view fname,
937 : bool useMolfile,
938 : int command_line_natoms) {
939 1264 : auto trajectory_fmt =FMTfromString(fmt);
940 1264 : if (useMolfile) {
941 : trajectory_fmt=trajfmt::molfile;
942 : }
943 901 : if (trajectory_fmt==trajfmt::error) {
944 : return "Error: wrong input format specified";
945 : }
946 :
947 2528 : parser = parserFromFMT(trajectory_fmt);
948 1264 : return parser->init(fmt,fname,command_line_natoms);
949 : }
950 :
951 78855 : std::optional<std::string> TrajectoryParser::readHeader(
952 : long long int &step,
953 : double &timeStep
954 : ) {
955 78855 : return parser->readHeader(step,timeStep);
956 : }
957 :
958 77616 : std::optional<std::string> TrajectoryParser::readAtoms(
959 : int stride,
960 : bool dont_read_pbc,
961 : bool debug_pd,
962 : int pd_start,
963 : int pd_nlocal,
964 : long long int &step,
965 : double* masses,
966 : double* charges,
967 : double* coordinates,
968 : double* cell
969 : ) {
970 77616 : return parser->readAtoms(
971 : stride,
972 : dont_read_pbc,
973 : debug_pd,
974 : pd_start,
975 : pd_nlocal,
976 : step,
977 : masses,
978 : charges,
979 : coordinates,
980 : cell
981 77616 : );
982 : }
983 :
984 139 : std::optional<std::string> TrajectoryParser::readFrame(
985 : int stride,
986 : bool dont_read_pbc,
987 : bool debug_pd,
988 : int pd_start,
989 : int pd_nlocal,
990 : long long int &step,
991 : double &timeStep,
992 : double * masses,
993 : double * charges,
994 : double * coordinates,
995 : double * cell
996 : ) {
997 139 : auto error = readHeader(step,timeStep);
998 139 : if (error) {
999 : return error;
1000 : }
1001 : ////////////////////////////////////////////////////////////////////////////
1002 : // here the first step settings
1003 : ////////////////////////////////////////////////////////////////////////////
1004 : //now the actual reading of the trajectory
1005 : return readAtoms(
1006 : stride,
1007 : dont_read_pbc,
1008 : debug_pd,
1009 : pd_start,
1010 : pd_nlocal,
1011 : step,
1012 : masses,
1013 : charges,
1014 : coordinates,
1015 : cell
1016 131 : );
1017 :
1018 : }
1019 :
1020 :
1021 0 : std::optional<std::string> TrajectoryParser::readFrame(
1022 : int stride,
1023 : bool dont_read_pbc,
1024 : bool debug_pd,
1025 : int pd_start,
1026 : int pd_nlocal,
1027 : long long int &step,
1028 : float &timeStep,
1029 : float * masses,
1030 : float * charges,
1031 : float * coordinates,
1032 : float * cell
1033 : ) {
1034 0 : auto error = readHeader(step,timeStep);
1035 0 : if (error) {
1036 : return error;
1037 : }
1038 : ////////////////////////////////////////////////////////////////////////////
1039 : // here the first step settings
1040 : ////////////////////////////////////////////////////////////////////////////
1041 : //now the actual reading of the trajectory
1042 : return readAtoms(
1043 : stride,
1044 : dont_read_pbc,
1045 : debug_pd,
1046 : pd_start,
1047 : pd_nlocal,
1048 : step,
1049 : masses,
1050 : charges,
1051 : coordinates,
1052 : cell
1053 0 : );
1054 : }
1055 :
1056 5 : std::optional<std::string> TrajectoryParser::readHeader(
1057 : long long int &step,
1058 : float &timeStep
1059 : ) {
1060 5 : return parser->readHeader(step,timeStep);
1061 : }
1062 :
1063 4 : std::optional<std::string> TrajectoryParser::readAtoms(
1064 : int stride,
1065 : bool dont_read_pbc,
1066 : bool debug_pd,
1067 : int pd_start,
1068 : int pd_nlocal,
1069 : long long int &step,
1070 : float* masses,
1071 : float* charges,
1072 : float* coordinates,
1073 : float* cell
1074 : ) {
1075 4 : return parser->readAtoms(stride,
1076 : dont_read_pbc,
1077 : debug_pd,
1078 : pd_start,
1079 : pd_nlocal,
1080 : step,
1081 : masses,
1082 : charges,
1083 : coordinates,
1084 4 : cell);
1085 : }
1086 :
1087 :
1088 : //some getters
1089 78791 : int TrajectoryParser::nOfAtoms() const {
1090 78791 : if (parser) {
1091 78791 : return parser->nOfAtoms();
1092 : }
1093 : return -1;
1094 : }
1095 :
1096 1 : std::optional<std::string> TrajectoryParser::rewind() {
1097 1 : return parser->rewind();
1098 : }
1099 :
1100 : } // namespace PLMD
|
