forked from pemsley/coot
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathextra-restraints.hh
269 lines (230 loc) · 8.83 KB
/
extra-restraints.hh
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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
/* ideal/extra-restraints.hh
*
* Copyright 2013 by Medical Research Council
* Author: Paul Emsley
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU General Public License and
* the GNU Lesser General Public License along with this program; if not,
* write to the Free Software Foundation, Inc., 51 Franklin Street,
* Fifth Floor, Boston, MA, 02110-1301, USA.
*/
#ifndef EXTRA_RESTRAINTS_HH
#define EXTRA_RESTRAINTS_HH
#include "coot-utils/coot-coord-utils.hh"
#include "parallel-planes.hh"
namespace coot {
class extra_restraints_t {
bool matches_bond_template_p(const std::vector<std::string> &words) const;
bool matches_angle_template_p(const std::vector<std::string> &words) const;
public:
class extra_bond_restraint_t {
public:
atom_spec_t atom_1;
atom_spec_t atom_2;
double bond_dist;
double esd;
extra_bond_restraint_t() {}
extra_bond_restraint_t(const atom_spec_t &a1, const atom_spec_t &a2, double d, double e) : atom_1(a1), atom_2(a2) {
bond_dist = d;
esd = e;
}
#ifdef SWIG
#else
bool operator==(const residue_spec_t &rs) const {
if (residue_spec_t(atom_1) == rs)
return true;
if (residue_spec_t(atom_2) == rs)
return true;
return false;
}
#endif
bool is_deviant(const double &real_dist, const double &n_sigma) const {
return (fabs(real_dist-bond_dist)/esd > n_sigma);
}
};
class bond_eraser {
public:
double n_sigma_lim;
std::map<std::pair<atom_spec_t, atom_spec_t>, double,
bool(*)(const std::pair<atom_spec_t, atom_spec_t> &,
const std::pair<atom_spec_t, atom_spec_t> &) > dist_map;
bond_eraser(const std::map<std::pair<atom_spec_t, atom_spec_t>, double,
bool(*)(const std::pair<atom_spec_t, atom_spec_t> &,
const std::pair<atom_spec_t, atom_spec_t> &) > &dist_map_in, double nsi) : dist_map(dist_map_in) {
n_sigma_lim = nsi;
}
bool operator()(const extra_bond_restraint_t &br) {
std::pair<atom_spec_t, atom_spec_t> p(br.atom_1, br.atom_2);
std::map<std::pair<atom_spec_t, atom_spec_t>, double>::const_iterator it =
dist_map.find(p);
if (it != dist_map.end()) {
bool v = (fabs(br.bond_dist - it->second)/br.esd >= n_sigma_lim);
std::cout << "comparing (" << br.bond_dist << " - " << it->second << ")/" << br.esd
<< " = " << fabs(br.bond_dist - it->second)/br.esd << " vs " << n_sigma_lim
<< " -> " << v << std::endl;
return (fabs(br.bond_dist - it->second)/br.esd >= n_sigma_lim);
} else {
std::cout << "not found bond restraint between " << br.atom_1 << " " << br.atom_2 << std::endl;
return false;
}
}
};
class extra_geman_mcclure_restraint_t : public extra_bond_restraint_t {
public:
extra_geman_mcclure_restraint_t() {}
extra_geman_mcclure_restraint_t(const atom_spec_t &a1, const atom_spec_t &a2, double d, double e) : extra_bond_restraint_t(a1, a2, d, e) {}
};
class extra_angle_restraint_t {
public:
atom_spec_t atom_1;
atom_spec_t atom_2;
atom_spec_t atom_3;
double angle;
double esd;
extra_angle_restraint_t(const atom_spec_t &a1, const atom_spec_t &a2,
const atom_spec_t &a3,
double angle_in, double esd_in) {
atom_1 = a1;
atom_2 = a2;
atom_3 = a3;
angle = angle_in;
esd = esd_in;
}
};
class extra_torsion_restraint_t {
public:
atom_spec_t atom_1;
atom_spec_t atom_2;
atom_spec_t atom_3;
atom_spec_t atom_4;
double torsion_angle;
double esd;
int period;
extra_torsion_restraint_t(const atom_spec_t &a1, const atom_spec_t &a2,
const atom_spec_t &a3, const atom_spec_t &a4,
double torsion_in, double esd_in, int period_in) {
atom_1 = a1;
atom_2 = a2;
atom_3 = a3;
atom_4 = a4;
torsion_angle = torsion_in;
esd = esd_in;
period = period_in;
}
};
class extra_start_pos_restraint_t {
public:
atom_spec_t atom_1;
double esd;
extra_start_pos_restraint_t(const atom_spec_t &a1, double e) {
atom_1 = a1;
esd = e;
}
};
class extra_target_position_restraint_t {
public:
atom_spec_t atom_spec;
clipper::Coord_orth pos;
double weight;
extra_target_position_restraint_t(const atom_spec_t &as, const clipper::Coord_orth &pos_in, double w) :
atom_spec(as), pos(pos_in), weight(w) {}
};
std::vector<extra_bond_restraint_t> bond_restraints;
std::vector<extra_angle_restraint_t> angle_restraints;
std::vector<extra_torsion_restraint_t> torsion_restraints;
std::vector<extra_geman_mcclure_restraint_t> geman_mcclure_restraints;
std::vector<extra_start_pos_restraint_t> start_pos_restraints;
std::vector<extra_target_position_restraint_t> target_position_restraints;
std::vector<parallel_planes_t> parallel_plane_restraints;
void read_refmac_extra_restraints(const std::string &file_name);
bool has_restraints() const {
if (bond_restraints.size() > 0)
return true;
if (angle_restraints.size() > 0)
return true;
else if (torsion_restraints.size() > 0)
return true;
else if (start_pos_restraints.size() > 0)
return true;
else if (parallel_plane_restraints.size() > 0)
return true;
else if (target_position_restraints.size() > 0)
return true;
else if (geman_mcclure_restraints.size() > 0)
return true;
else
return false;
}
void clear() {
bond_restraints.clear();
angle_restraints.clear();
torsion_restraints.clear();
start_pos_restraints.clear();
parallel_plane_restraints.clear();
geman_mcclure_restraints.clear();
}
void add_restraints(const extra_restraints_t &r) {
for (unsigned int i=0; i<r.bond_restraints.size(); i++)
bond_restraints.push_back(r.bond_restraints[i]);
for (unsigned int i=0; i<r.angle_restraints.size(); i++)
angle_restraints.push_back(r.angle_restraints[i]);
for (unsigned int i=0; i<r.torsion_restraints.size(); i++)
torsion_restraints.push_back(r.torsion_restraints[i]);
for (unsigned int i=0; i<r.start_pos_restraints.size(); i++)
start_pos_restraints.push_back(r.start_pos_restraints[i]);
for (unsigned int i=0; i<r.parallel_plane_restraints.size(); i++)
parallel_plane_restraints.push_back(r.parallel_plane_restraints[i]);
// 20231003-PE GM was not in this list. Why not?
for (unsigned int i=0; i<r.geman_mcclure_restraints.size(); i++)
geman_mcclure_restraints.push_back(r.geman_mcclure_restraints[i]);
}
void delete_restraints_for_residue(const residue_spec_t &rs);
// updates restraint on atom if it can, else adds
void add_start_pos_restraint(const atom_spec_t &atom_1_in, double esd_in);
// We want to interpolate proSMART restraints from start to final model.
// We have proSMART restraints for both models.
//
// Let's return a list of extra bond restraint indices that are
// between the bond restraints of this (presumably start) and
// final.
//
std::vector<std::pair<unsigned int, unsigned int> >
find_pair_indices(const extra_restraints_t &final) const;
// n_path_points is the number of points along the trajectory.
// It should be at least 2.
//
void
write_interpolated_restraints(const extra_restraints_t &final,
unsigned int n_path_points,
std::string file_name_stub) const;
std::map<mmdb::Atom *, clipper::Coord_orth>
position_point_map(mmdb::Manager *mol_running, mmdb::Manager *mol_ref) const;
void
write_interpolated_models_and_restraints(const extra_restraints_t &final,
mmdb::Manager *mol_1, // corresponds to this
mmdb::Manager *mol_2, // corresponds to final
unsigned int n_path_points,
std::string file_name_stub) const;
void write_interpolated_restraints(std::ofstream &f,
const std::vector<extra_bond_restraint_t> &final_bonds_restraints,
double frac,
unsigned int idx_1,
unsigned int idx2) const;
void write_interpolated_models(mmdb::Manager *mol_running,
const std::map<mmdb::Atom *, clipper::Coord_orth> &matching_atoms_1,
const std::map<mmdb::Atom *, clipper::Coord_orth> &matching_atoms_2,
unsigned int n_path_points,
std::string file_name_stub) const;
};
}
#endif // EXTRA_RESTRAINTS_HH