1 #include <cmath>
2
3 #include "colvarmodule.h"
4 #include "colvarvalue.h"
5 #include "colvarparse.h"
6 #include "colvar.h"
7 #include "colvarcomp.h"
8
9
10
orientation(std::string const & conf)11 colvar::orientation::orientation (std::string const &conf)
12 : cvc (conf)
13 {
14 function_type = "orientation";
15 parse_group (conf, "atoms", atoms);
16 atom_groups.push_back (&atoms);
17 x.type (colvarvalue::type_quaternion);
18
19 ref_pos.reserve (atoms.size());
20
21 if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
22 cvm::log ("Using reference positions from input file.\n");
23 if (ref_pos.size() != atoms.size()) {
24 cvm::fatal_error ("Error: reference positions do not "
25 "match the number of requested atoms.\n");
26 }
27 }
28
29 {
30 std::string file_name;
31 if (get_keyval (conf, "refPositionsFile", file_name)) {
32
33 std::string file_col;
34 double file_col_value;
35 if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) {
36 // use PDB flags if column is provided
37 bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0);
38 if (found && !file_col_value)
39 cvm::fatal_error ("Error: refPositionsColValue, "
40 "if provided, must be non-zero.\n");
41 } else {
42 // if not, use atom indices
43 atoms.create_sorted_ids();
44 }
45 ref_pos.resize (atoms.size());
46 cvm::load_coords (file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value);
47 }
48 }
49
50 if (!ref_pos.size()) {
51 cvm::fatal_error ("Error: must define a set of "
52 "reference coordinates.\n");
53 }
54
55
56 cvm::log ("Centering the reference coordinates: it is "
57 "assumed that each atom is the closest "
58 "periodic image to the center of geometry.\n");
59 cvm::rvector cog (0.0, 0.0, 0.0);
60 for (size_t i = 0; i < ref_pos.size(); i++) {
61 cog += ref_pos[i];
62 }
63 cog /= cvm::real (ref_pos.size());
64 for (size_t i = 0; i < ref_pos.size(); i++) {
65 ref_pos[i] -= cog;
66 }
67
68 get_keyval (conf, "closestToQuaternion", ref_quat, cvm::quaternion (1.0, 0.0, 0.0, 0.0));
69
70 // initialize rot member data
71 if (!atoms.noforce) {
72 rot.request_group2_gradients (atoms.size());
73 }
74
75 }
76
77
orientation()78 colvar::orientation::orientation()
79 : cvc ()
80 {
81 function_type = "orientation";
82 x.type (colvarvalue::type_quaternion);
83 }
84
85
calc_value()86 void colvar::orientation::calc_value()
87 {
88 // atoms.reset_atoms_data();
89 // atoms.read_positions();
90
91 atoms_cog = atoms.center_of_geometry();
92
93 rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
94
95 if ((rot.q).inner (ref_quat) >= 0.0) {
96 x.quaternion_value = rot.q;
97 } else {
98 x.quaternion_value = -1.0 * rot.q;
99 }
100 }
101
102
calc_gradients()103 void colvar::orientation::calc_gradients()
104 {
105 // gradients have already been calculated and stored within the
106 // member object "rot"; we're not using the "grad" member of each
107 // atom object, because it only can represent the gradient of a
108 // scalar colvar
109 }
110
111
apply_force(colvarvalue const & force)112 void colvar::orientation::apply_force (colvarvalue const &force)
113 {
114 cvm::quaternion const &FQ = force.quaternion_value;
115
116 if (!atoms.noforce) {
117 for (size_t ia = 0; ia < atoms.size(); ia++) {
118 for (size_t i = 0; i < 4; i++) {
119 atoms[ia].apply_force (FQ[i] * rot.dQ0_2[ia][i]);
120 }
121 }
122 }
123 }
124
125
126
orientation_angle(std::string const & conf)127 colvar::orientation_angle::orientation_angle (std::string const &conf)
128 : orientation (conf)
129 {
130 function_type = "orientation_angle";
131 x.type (colvarvalue::type_scalar);
132 }
133
134
orientation_angle()135 colvar::orientation_angle::orientation_angle()
136 : orientation()
137 {
138 function_type = "orientation_angle";
139 x.type (colvarvalue::type_scalar);
140 }
141
142
calc_value()143 void colvar::orientation_angle::calc_value()
144 {
145 // atoms.reset_atoms_data();
146 // atoms.read_positions();
147
148 atoms_cog = atoms.center_of_geometry();
149
150 rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
151
152 if ((rot.q).q0 >= 0.0) {
153 x.real_value = (180.0/PI) * 2.0 * std::acos ((rot.q).q0);
154 } else {
155 x.real_value = (180.0/PI) * 2.0 * std::acos (-1.0 * (rot.q).q0);
156 }
157 }
158
159
calc_gradients()160 void colvar::orientation_angle::calc_gradients()
161 {
162 cvm::real const dxdq0 =
163 ( ((rot.q).q0 * (rot.q).q0 < 1.0) ?
164 ((180.0 / PI) * (-2.0) / std::sqrt (1.0 - ((rot.q).q0 * (rot.q).q0))) :
165 0.0 );
166
167 for (size_t ia = 0; ia < atoms.size(); ia++) {
168 atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
169 }
170 }
171
172
apply_force(colvarvalue const & force)173 void colvar::orientation_angle::apply_force (colvarvalue const &force)
174 {
175 cvm::real const &fw = force.real_value;
176
177 if (!atoms.noforce) {
178 atoms.apply_colvar_force (fw);
179 }
180 }
181
182
tilt(std::string const & conf)183 colvar::tilt::tilt (std::string const &conf)
184 : orientation (conf)
185 {
186 function_type = "tilt";
187
188 get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0));
189
190 if (axis.norm2() != 1.0) {
191 axis /= axis.norm();
192 cvm::log ("Normalizing rotation axis to "+cvm::to_str (axis)+".\n");
193 }
194
195 x.type (colvarvalue::type_scalar);
196 }
197
198
tilt()199 colvar::tilt::tilt()
200 : orientation()
201 {
202 function_type = "tilt";
203 x.type (colvarvalue::type_scalar);
204 }
205
206
calc_value()207 void colvar::tilt::calc_value()
208 {
209 // atoms.reset_atoms_data();
210 // atoms.read_positions();
211
212 atoms_cog = atoms.center_of_geometry();
213
214 rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
215
216 x.real_value = rot.cos_theta (axis);
217 }
218
219
calc_gradients()220 void colvar::tilt::calc_gradients()
221 {
222 cvm::quaternion const dxdq = rot.dcos_theta_dq (axis);
223
224 for (size_t ia = 0; ia < atoms.size(); ia++) {
225 atoms[ia].grad = cvm::rvector (0.0, 0.0, 0.0);
226 for (size_t iq = 0; iq < 4; iq++) {
227 atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
228 }
229 }
230
231 if (b_debug_gradients) {
232 cvm::log ("Debugging tilt component gradients:\n");
233 debug_gradients (atoms);
234 }
235 }
236
237
apply_force(colvarvalue const & force)238 void colvar::tilt::apply_force (colvarvalue const &force)
239 {
240 cvm::real const &fw = force.real_value;
241
242 if (!atoms.noforce) {
243 atoms.apply_colvar_force (fw);
244 }
245 }
246
247
248
spin_angle(std::string const & conf)249 colvar::spin_angle::spin_angle (std::string const &conf)
250 : orientation (conf)
251 {
252 function_type = "spin_angle";
253
254 get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0));
255
256 if (axis.norm2() != 1.0) {
257 axis /= axis.norm();
258 cvm::log ("Normalizing rotation axis to "+cvm::to_str (axis)+".\n");
259 }
260
261 period = 360.0;
262 b_periodic = true;
263 x.type (colvarvalue::type_scalar);
264 }
265
266
spin_angle()267 colvar::spin_angle::spin_angle()
268 : orientation()
269 {
270 function_type = "spin_angle";
271 period = 360.0;
272 b_periodic = true;
273 x.type (colvarvalue::type_scalar);
274 }
275
276
calc_value()277 void colvar::spin_angle::calc_value()
278 {
279 // atoms.reset_atoms_data();
280 // atoms.read_positions();
281
282 atoms_cog = atoms.center_of_geometry();
283
284 rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
285
286 x.real_value = rot.spin_angle (axis);
287 this->wrap (x);
288 }
289
290
calc_gradients()291 void colvar::spin_angle::calc_gradients()
292 {
293 cvm::quaternion const dxdq = rot.dspin_angle_dq (axis);
294
295 for (size_t ia = 0; ia < atoms.size(); ia++) {
296 atoms[ia].grad = cvm::rvector (0.0, 0.0, 0.0);
297 for (size_t iq = 0; iq < 4; iq++) {
298 atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
299 }
300 }
301 }
302
303
apply_force(colvarvalue const & force)304 void colvar::spin_angle::apply_force (colvarvalue const &force)
305 {
306 cvm::real const &fw = force.real_value;
307
308 if (!atoms.noforce) {
309 atoms.apply_colvar_force (fw);
310 }
311 }
312