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