1 //
2 // Copyright (C) 2003-2016 Sereina Riniker, Paolo Tosco
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10
11 #include <RDGeneral/Invariant.h>
12 #include <RDGeneral/BadFileException.h>
13 #include <GraphMol/ROMol.h>
14 #include <GraphMol/Conformer.h>
15 #include <fstream>
16 #include <sstream>
17 #include <set>
18 #include "Trajectory.h"
19
20 namespace RDKit {
21
getPoint2D(unsigned int pointNum) const22 RDGeom::Point2D Snapshot::getPoint2D(unsigned int pointNum) const {
23 PRECONDITION(d_pos, "d_pos must not be NULL");
24 PRECONDITION(d_trajectory, "d_trajectory must not be NULL");
25 PRECONDITION(d_trajectory->dimension() == 2, "d_dimension must be == 2");
26 PRECONDITION(d_trajectory->numPoints(), "d_numPoints must be > 0");
27 URANGE_CHECK(pointNum, d_trajectory->numPoints());
28 unsigned int i = pointNum * d_trajectory->dimension();
29 return RDGeom::Point2D(d_pos[i], d_pos[i + 1]);
30 }
31
getPoint3D(unsigned int pointNum) const32 RDGeom::Point3D Snapshot::getPoint3D(unsigned int pointNum) const {
33 PRECONDITION(d_pos, "d_pos must not be NULL");
34 PRECONDITION(d_trajectory, "d_trajectory must not be NULL");
35 PRECONDITION(d_trajectory->dimension() >= 2, "d_dimension must be >= 2");
36 PRECONDITION(d_trajectory->numPoints(), "d_numPoints must be > 0");
37 URANGE_CHECK(pointNum, d_trajectory->numPoints());
38 unsigned int i = pointNum * d_trajectory->dimension();
39 return (
40 RDGeom::Point3D(d_pos[i], d_pos[i + 1],
41 (d_trajectory->dimension() == 3) ? d_pos[i + 2] : 0.0));
42 }
43
Trajectory(unsigned int dimension,unsigned int numPoints,SnapshotVect * snapshotVect)44 Trajectory::Trajectory(unsigned int dimension, unsigned int numPoints,
45 SnapshotVect *snapshotVect)
46 : d_dimension(dimension), d_numPoints(numPoints) {
47 if (!snapshotVect) {
48 snapshotVect = new SnapshotVect;
49 }
50 d_snapshotVect.reset(snapshotVect);
51 for (auto &vectIt : *d_snapshotVect) {
52 vectIt.d_trajectory = this;
53 }
54 }
55
Trajectory(const Trajectory & other)56 Trajectory::Trajectory(const Trajectory &other)
57 : d_dimension(other.d_dimension),
58 d_numPoints(other.d_numPoints),
59 d_snapshotVect(new SnapshotVect) {
60 for (SnapshotVect::const_iterator vectIt = other.d_snapshotVect->begin();
61 vectIt != other.d_snapshotVect->end(); ++vectIt) {
62 addSnapshot(*vectIt);
63 }
64 }
65
addSnapshot(const Snapshot & s)66 unsigned int Trajectory::addSnapshot(const Snapshot &s) {
67 return insertSnapshot(d_snapshotVect->size(), s);
68 }
69
getSnapshot(unsigned int snapshotNum) const70 const Snapshot &Trajectory::getSnapshot(unsigned int snapshotNum) const {
71 URANGE_CHECK(snapshotNum, d_snapshotVect->size());
72 return (*d_snapshotVect)[snapshotNum];
73 }
74
insertSnapshot(unsigned int snapshotNum,Snapshot s)75 unsigned int Trajectory::insertSnapshot(unsigned int snapshotNum, Snapshot s) {
76 URANGE_CHECK(snapshotNum, d_snapshotVect->size() + 1);
77 s.d_trajectory = this;
78 return (d_snapshotVect->insert(d_snapshotVect->begin() + snapshotNum, s) -
79 d_snapshotVect->begin());
80 }
81
removeSnapshot(unsigned int snapshotNum)82 unsigned int Trajectory::removeSnapshot(unsigned int snapshotNum) {
83 URANGE_CHECK(snapshotNum, d_snapshotVect->size());
84 return (d_snapshotVect->erase(d_snapshotVect->begin() + snapshotNum) -
85 d_snapshotVect->begin());
86 }
87
addConformersToMol(ROMol & mol,int from,int to)88 unsigned int Trajectory::addConformersToMol(ROMol &mol, int from, int to) {
89 PRECONDITION(d_numPoints == mol.getNumAtoms(),
90 "Number of atom mismatch between ROMol and Trajectory");
91 PRECONDITION(from < static_cast<int>(size()), "from must be < size()");
92 PRECONDITION(to < static_cast<int>(size()), "to must be < size()");
93 if (from < 0) {
94 from = 0;
95 }
96 if (to < 0) {
97 to = size() - 1;
98 }
99 PRECONDITION(!size() || (from <= to), "from must be <= to");
100 int n;
101 unsigned int nConf;
102 for (n = from, nConf = 0; size() && (n <= to); ++n, ++nConf) {
103 auto *conf = new Conformer(mol.getNumAtoms());
104 for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
105 conf->setAtomPos(i, getSnapshot(n).getPoint3D(i));
106 }
107 mol.addConformer(conf, true);
108 }
109 return nConf;
110 }
111
readAmberTrajectory(const std::string & fName,Trajectory & traj)112 unsigned int readAmberTrajectory(const std::string &fName, Trajectory &traj) {
113 PRECONDITION(traj.dimension() == 3,
114 "The trajectory must have dimension == 3");
115 std::ifstream inStream(fName.c_str());
116 if (!inStream || inStream.bad()) {
117 std::stringstream ss;
118 ss << "Bad input file: " << fName;
119 throw RDKit::BadFileException(ss.str());
120 }
121 std::string tempStr;
122 // title
123 std::getline(inStream, tempStr);
124 // read coordinates
125 unsigned int nCoords = traj.numPoints() * 3;
126 unsigned int nSnapshots = 0;
127 while (inStream.good() && !inStream.eof()) {
128 boost::shared_array<double> c(new double[nCoords]());
129 unsigned int i = 0;
130 while (i < nCoords) {
131 if (!(inStream >> c[i])) {
132 if (!inStream.eof()) {
133 std::stringstream ss;
134 ss << "Error while reading file: " << fName;
135 throw ValueErrorException(ss.str());
136 } else if (i && (i < (nCoords - 1))) {
137 std::stringstream ss;
138 ss << "Premature end of file: " << fName;
139 throw ValueErrorException(ss.str());
140 }
141 break;
142 }
143 ++i;
144 }
145 if (!inStream.eof()) {
146 traj.addSnapshot(Snapshot(c));
147 ++nSnapshots;
148 }
149 }
150 return nSnapshots;
151 }
152
readGromosTrajectory(const std::string & fName,Trajectory & traj)153 unsigned int readGromosTrajectory(const std::string &fName, Trajectory &traj) {
154 PRECONDITION(traj.dimension() == 3,
155 "The trajectory must have dimension == 3");
156 std::ifstream inStream(fName.c_str());
157 if (!inStream || inStream.bad()) {
158 std::stringstream ss;
159 ss << "Bad input file: " << fName;
160 throw RDKit::BadFileException(ss.str());
161 }
162 std::string tempStr;
163 unsigned int nCoords = traj.numPoints() * 3;
164 unsigned int nSnapshots = 0;
165 const static char *ignoredKeywordArray[] = {
166 "TITLE", "TIMESTEP", "VELOCITYRED", "VELOCITY", "GENBOX", "BOX"};
167 std::set<std::string> ignoredKeywordSet;
168 for (auto &kw : ignoredKeywordArray) {
169 ignoredKeywordSet.insert(std::string(kw));
170 }
171 while (inStream.good() && !inStream.eof()) {
172 std::getline(inStream, tempStr);
173 if (inStream.bad() || inStream.eof()) {
174 continue;
175 }
176 if (ignoredKeywordSet.find(tempStr) != ignoredKeywordSet.end()) {
177 // ignored block
178 while (inStream.good() && !inStream.eof() && (tempStr != "END")) {
179 std::getline(inStream, tempStr);
180 }
181 } else if ((tempStr == "POSITIONRED") || (tempStr == "POSITION")) {
182 // these are the positions
183 boost::shared_array<double> c(new double[nCoords]());
184 unsigned int j = 0;
185 for (unsigned int i = 0; i < traj.numPoints();) {
186 std::getline(inStream, tempStr);
187 if (inStream.bad() || inStream.eof() || (tempStr == "END")) {
188 throw ValueErrorException("Wrong number of coordinates");
189 }
190 // ignore comments
191 if (tempStr.find("#") != std::string::npos) {
192 continue;
193 }
194 std::stringstream ls(tempStr);
195 double x, y, z;
196 if (!(ls >> x >> y >> z)) {
197 throw ValueErrorException("Error while reading file");
198 }
199 // store the coordinates (convert to Angstrom!)
200 c[j++] = x * 10.0;
201 c[j++] = y * 10.0;
202 c[j++] = z * 10.0;
203 ++i;
204 }
205 std::getline(inStream, tempStr); // the END line
206 if (inStream.bad() || inStream.eof() || (tempStr != "END")) {
207 throw ValueErrorException("Wrong number of coordinates");
208 }
209 traj.addSnapshot(Snapshot(c));
210 ++nSnapshots;
211 } else {
212 std::string supportedBlocks("POSITIONRED, POSITION");
213 for (const auto &it : ignoredKeywordSet) {
214 supportedBlocks += ", " + it;
215 }
216 throw ValueErrorException("Unsupported block: " + tempStr +
217 ". Supported blocks are " + supportedBlocks);
218 }
219 } // read file
220 if (inStream.bad()) {
221 std::stringstream ss;
222 ss << "Bad input file: " << fName;
223 throw RDKit::BadFileException(ss.str());
224 }
225 return nSnapshots;
226 }
227 }
228