1 // -*- c++ -*-
2
3 // This file is part of the Collective Variables module (Colvars).
4 // The original version of Colvars and its updates are located at:
5 // https://github.com/Colvars/colvars
6 // Please update all Colvars source files before making any changes.
7 // If you wish to distribute your changes, please submit them to the
8 // Colvars repository at GitHub.
9
10 #include "colvarmodule.h"
11 #include "colvar.h"
12 #include "colvarcomp.h"
13
14
angle(std::string const & conf)15 colvar::angle::angle(std::string const &conf)
16 : cvc(conf)
17 {
18 function_type = "angle";
19 init_as_angle();
20
21 provide(f_cvc_inv_gradient);
22 provide(f_cvc_Jacobian);
23 enable(f_cvc_com_based);
24
25 group1 = parse_group(conf, "group1");
26 group2 = parse_group(conf, "group2");
27 group3 = parse_group(conf, "group3");
28
29 init_total_force_params(conf);
30 }
31
32
angle(cvm::atom const & a1,cvm::atom const & a2,cvm::atom const & a3)33 colvar::angle::angle(cvm::atom const &a1,
34 cvm::atom const &a2,
35 cvm::atom const &a3)
36 {
37 function_type = "angle";
38 init_as_angle();
39
40 provide(f_cvc_inv_gradient);
41 provide(f_cvc_Jacobian);
42 enable(f_cvc_com_based);
43
44 group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
45 group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
46 group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
47 register_atom_group(group1);
48 register_atom_group(group2);
49 register_atom_group(group3);
50 }
51
52
calc_value()53 void colvar::angle::calc_value()
54 {
55 cvm::atom_pos const g1_pos = group1->center_of_mass();
56 cvm::atom_pos const g2_pos = group2->center_of_mass();
57 cvm::atom_pos const g3_pos = group3->center_of_mass();
58
59 r21 = is_enabled(f_cvc_pbc_minimum_image) ?
60 cvm::position_distance(g2_pos, g1_pos) :
61 g1_pos - g2_pos;
62 r21l = r21.norm();
63 r23 = is_enabled(f_cvc_pbc_minimum_image) ?
64 cvm::position_distance(g2_pos, g3_pos) :
65 g3_pos - g2_pos;
66 r23l = r23.norm();
67
68 cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
69
70 x.real_value = (180.0/PI) * cvm::acos(cos_theta);
71 }
72
73
calc_gradients()74 void colvar::angle::calc_gradients()
75 {
76 cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
77 cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta);
78
79 dxdr1 = (180.0/PI) * dxdcos *
80 (1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l );
81
82 dxdr3 = (180.0/PI) * dxdcos *
83 (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
84
85 group1->set_weighted_gradient(dxdr1);
86 group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0));
87 group3->set_weighted_gradient(dxdr3);
88 }
89
90
calc_force_invgrads()91 void colvar::angle::calc_force_invgrads()
92 {
93 // This uses a force measurement on groups 1 and 3 only
94 // to keep in line with the implicit variable change used to
95 // evaluate the Jacobian term (essentially polar coordinates
96 // centered on group2, which means group2 is kept fixed
97 // when propagating changes in the angle)
98
99 if (is_enabled(f_cvc_one_site_total_force)) {
100 group1->read_total_forces();
101 cvm::real norm_fact = 1.0 / dxdr1.norm2();
102 ft.real_value = norm_fact * dxdr1 * group1->total_force();
103 } else {
104 group1->read_total_forces();
105 group3->read_total_forces();
106 cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2());
107 ft.real_value = norm_fact * ( dxdr1 * group1->total_force()
108 + dxdr3 * group3->total_force());
109 }
110 return;
111 }
112
113
calc_Jacobian_derivative()114 void colvar::angle::calc_Jacobian_derivative()
115 {
116 // det(J) = (2 pi) r^2 * sin(theta)
117 // hence Jd = cot(theta)
118 const cvm::real theta = x.real_value * PI / 180.0;
119 jd = PI / 180.0 * (theta != 0.0 ? cvm::cos(theta) / cvm::sin(theta) : 0.0);
120 }
121
122
apply_force(colvarvalue const & force)123 void colvar::angle::apply_force(colvarvalue const &force)
124 {
125 if (!group1->noforce)
126 group1->apply_colvar_force(force.real_value);
127
128 if (!group2->noforce)
129 group2->apply_colvar_force(force.real_value);
130
131 if (!group3->noforce)
132 group3->apply_colvar_force(force.real_value);
133 }
134
135
simple_scalar_dist_functions(angle)136 simple_scalar_dist_functions(angle)
137
138
139
140 colvar::dipole_angle::dipole_angle(std::string const &conf)
141 : cvc(conf)
142 {
143 function_type = "dipole_angle";
144 init_as_angle();
145
146 group1 = parse_group(conf, "group1");
147 group2 = parse_group(conf, "group2");
148 group3 = parse_group(conf, "group3");
149
150 init_total_force_params(conf);
151 }
152
153
dipole_angle(cvm::atom const & a1,cvm::atom const & a2,cvm::atom const & a3)154 colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
155 cvm::atom const &a2,
156 cvm::atom const &a3)
157 {
158 function_type = "dipole_angle";
159
160 group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
161 group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
162 group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
163 register_atom_group(group1);
164 register_atom_group(group2);
165 register_atom_group(group3);
166
167 x.type(colvarvalue::type_scalar);
168 }
169
170
dipole_angle()171 colvar::dipole_angle::dipole_angle()
172 {
173 function_type = "dipole_angle";
174 init_as_angle();
175 x.type(colvarvalue::type_scalar);
176 }
177
178
calc_value()179 void colvar::dipole_angle::calc_value()
180 {
181 cvm::atom_pos const g1_pos = group1->center_of_mass();
182 cvm::atom_pos const g2_pos = group2->center_of_mass();
183 cvm::atom_pos const g3_pos = group3->center_of_mass();
184
185 group1->calc_dipole(g1_pos);
186
187 r21 = group1->dipole();
188 r21l = r21.norm();
189 r23 = is_enabled(f_cvc_pbc_minimum_image) ?
190 cvm::position_distance(g2_pos, g3_pos) :
191 g3_pos - g2_pos;
192 r23l = r23.norm();
193
194 cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
195
196 x.real_value = (180.0/PI) * cvm::acos(cos_theta);
197 }
198
199 //to be implemented
200 //void colvar::dipole_angle::calc_force_invgrads(){}
201 //void colvar::dipole_angle::calc_Jacobian_derivative(){}
202
calc_gradients()203 void colvar::dipole_angle::calc_gradients()
204 {
205 cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
206 cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta);
207
208 dxdr1 = (180.0/PI) * dxdcos *
209 (1.0/r21l)* (r23/r23l + (-1.0) * cos_theta * r21/r21l );
210
211 dxdr3 = (180.0/PI) * dxdcos *
212 (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
213
214 //this auxiliar variables are to avoid numerical errors inside "for"
215 double aux1 = group1->total_charge/group1->total_mass;
216 // double aux2 = group2->total_charge/group2->total_mass;
217 // double aux3 = group3->total_charge/group3->total_mass;
218
219 size_t i;
220 for (i = 0; i < group1->size(); i++) {
221 (*group1)[i].grad =((*group1)[i].charge + (-1)* (*group1)[i].mass * aux1) * (dxdr1);
222 }
223
224 for (i = 0; i < group2->size(); i++) {
225 (*group2)[i].grad = ((*group2)[i].mass/group2->total_mass)* dxdr3 * (-1.0);
226 }
227
228 for (i = 0; i < group3->size(); i++) {
229 (*group3)[i].grad =((*group3)[i].mass/group3->total_mass) * (dxdr3);
230 }
231 }
232
233
apply_force(colvarvalue const & force)234 void colvar::dipole_angle::apply_force(colvarvalue const &force)
235 {
236 if (!group1->noforce)
237 group1->apply_colvar_force(force.real_value);
238
239 if (!group2->noforce)
240 group2->apply_colvar_force(force.real_value);
241
242 if (!group3->noforce)
243 group3->apply_colvar_force(force.real_value);
244 }
245
246
simple_scalar_dist_functions(dipole_angle)247 simple_scalar_dist_functions(dipole_angle)
248
249
250
251 colvar::dihedral::dihedral(std::string const &conf)
252 : cvc(conf)
253 {
254 function_type = "dihedral";
255 period = 360.0;
256 enable(f_cvc_periodic);
257 provide(f_cvc_inv_gradient);
258 provide(f_cvc_Jacobian);
259 enable(f_cvc_com_based);
260
261 group1 = parse_group(conf, "group1");
262 group2 = parse_group(conf, "group2");
263 group3 = parse_group(conf, "group3");
264 group4 = parse_group(conf, "group4");
265
266 init_total_force_params(conf);
267
268 x.type(colvarvalue::type_scalar);
269 }
270
271
dihedral(cvm::atom const & a1,cvm::atom const & a2,cvm::atom const & a3,cvm::atom const & a4)272 colvar::dihedral::dihedral(cvm::atom const &a1,
273 cvm::atom const &a2,
274 cvm::atom const &a3,
275 cvm::atom const &a4)
276 {
277 if (cvm::debug())
278 cvm::log("Initializing dihedral object from atom groups.\n");
279
280 function_type = "dihedral";
281 period = 360.0;
282 enable(f_cvc_periodic);
283 provide(f_cvc_inv_gradient);
284 provide(f_cvc_Jacobian);
285 enable(f_cvc_com_based);
286
287 b_1site_force = false;
288
289 group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
290 group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
291 group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
292 group4 = new cvm::atom_group(std::vector<cvm::atom>(1, a4));
293 register_atom_group(group1);
294 register_atom_group(group2);
295 register_atom_group(group3);
296 register_atom_group(group4);
297
298 x.type(colvarvalue::type_scalar);
299
300 if (cvm::debug())
301 cvm::log("Done initializing dihedral object from atom groups.\n");
302 }
303
304
dihedral()305 colvar::dihedral::dihedral()
306 {
307 function_type = "dihedral";
308 period = 360.0;
309 enable(f_cvc_periodic);
310 provide(f_cvc_inv_gradient);
311 provide(f_cvc_Jacobian);
312 x.type(colvarvalue::type_scalar);
313 }
314
315
calc_value()316 void colvar::dihedral::calc_value()
317 {
318 cvm::atom_pos const g1_pos = group1->center_of_mass();
319 cvm::atom_pos const g2_pos = group2->center_of_mass();
320 cvm::atom_pos const g3_pos = group3->center_of_mass();
321 cvm::atom_pos const g4_pos = group4->center_of_mass();
322
323 // Usual sign convention: r12 = r2 - r1
324 r12 = is_enabled(f_cvc_pbc_minimum_image) ?
325 cvm::position_distance(g1_pos, g2_pos) :
326 g2_pos - g1_pos;
327 r23 = is_enabled(f_cvc_pbc_minimum_image) ?
328 cvm::position_distance(g2_pos, g3_pos) :
329 g3_pos - g2_pos;
330 r34 = is_enabled(f_cvc_pbc_minimum_image) ?
331 cvm::position_distance(g3_pos, g4_pos) :
332 g4_pos - g3_pos;
333
334 cvm::rvector const n1 = cvm::rvector::outer(r12, r23);
335 cvm::rvector const n2 = cvm::rvector::outer(r23, r34);
336
337 cvm::real const cos_phi = n1 * n2;
338 cvm::real const sin_phi = n1 * r34 * r23.norm();
339
340 x.real_value = (180.0/PI) * cvm::atan2(sin_phi, cos_phi);
341 this->wrap(x);
342 }
343
344
calc_gradients()345 void colvar::dihedral::calc_gradients()
346 {
347 cvm::rvector A = cvm::rvector::outer(r12, r23);
348 cvm::real rA = A.norm();
349 cvm::rvector B = cvm::rvector::outer(r23, r34);
350 cvm::real rB = B.norm();
351 cvm::rvector C = cvm::rvector::outer(r23, A);
352 cvm::real rC = C.norm();
353
354 cvm::real const cos_phi = (A*B)/(rA*rB);
355 cvm::real const sin_phi = (C*B)/(rC*rB);
356
357 cvm::rvector f1, f2, f3;
358
359 rB = 1.0/rB;
360 B *= rB;
361
362 if (cvm::fabs(sin_phi) > 0.1) {
363 rA = 1.0/rA;
364 A *= rA;
365 cvm::rvector const dcosdA = rA*(cos_phi*A-B);
366 cvm::rvector const dcosdB = rB*(cos_phi*B-A);
367 rA = 1.0;
368
369 cvm::real const K = (1.0/sin_phi) * (180.0/PI);
370
371 f1 = K * cvm::rvector::outer(r23, dcosdA);
372 f3 = K * cvm::rvector::outer(dcosdB, r23);
373 f2 = K * (cvm::rvector::outer(dcosdA, r12)
374 + cvm::rvector::outer(r34, dcosdB));
375 }
376 else {
377 rC = 1.0/rC;
378 C *= rC;
379 cvm::rvector const dsindC = rC*(sin_phi*C-B);
380 cvm::rvector const dsindB = rB*(sin_phi*B-C);
381 rC = 1.0;
382
383 cvm::real const K = (-1.0/cos_phi) * (180.0/PI);
384
385 f1.x = K*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
386 - r23.x*r23.y*dsindC.y
387 - r23.x*r23.z*dsindC.z);
388 f1.y = K*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
389 - r23.y*r23.z*dsindC.z
390 - r23.y*r23.x*dsindC.x);
391 f1.z = K*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
392 - r23.z*r23.x*dsindC.x
393 - r23.z*r23.y*dsindC.y);
394
395 f3 = cvm::rvector::outer(dsindB, r23);
396 f3 *= K;
397
398 f2.x = K*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
399 +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
400 +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
401 +dsindB.z*r34.y - dsindB.y*r34.z);
402 f2.y = K*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
403 +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
404 +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
405 +dsindB.x*r34.z - dsindB.z*r34.x);
406 f2.z = K*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
407 +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
408 +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
409 +dsindB.y*r34.x - dsindB.x*r34.y);
410 }
411
412 group1->set_weighted_gradient(-f1);
413 group2->set_weighted_gradient(-f2 + f1);
414 group3->set_weighted_gradient(-f3 + f2);
415 group4->set_weighted_gradient(f3);
416 }
417
418
calc_force_invgrads()419 void colvar::dihedral::calc_force_invgrads()
420 {
421 cvm::rvector const u12 = r12.unit();
422 cvm::rvector const u23 = r23.unit();
423 cvm::rvector const u34 = r34.unit();
424
425 cvm::real const d12 = r12.norm();
426 cvm::real const d34 = r34.norm();
427
428 cvm::rvector const cross1 = (cvm::rvector::outer(u23, u12)).unit();
429 cvm::rvector const cross4 = (cvm::rvector::outer(u23, u34)).unit();
430
431 cvm::real const dot1 = u23 * u12;
432 cvm::real const dot4 = u23 * u34;
433
434 cvm::real const fact1 = d12 * cvm::sqrt(1.0 - dot1 * dot1);
435 cvm::real const fact4 = d34 * cvm::sqrt(1.0 - dot4 * dot4);
436
437 group1->read_total_forces();
438 if (is_enabled(f_cvc_one_site_total_force)) {
439 // This is only measuring the force on group 1
440 ft.real_value = PI/180.0 * fact1 * (cross1 * group1->total_force());
441 } else {
442 // Default case: use groups 1 and 4
443 group4->read_total_forces();
444 ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1->total_force())
445 + fact4 * (cross4 * group4->total_force()));
446 }
447 }
448
449
calc_Jacobian_derivative()450 void colvar::dihedral::calc_Jacobian_derivative()
451 {
452 // With this choice of inverse gradient ("internal coordinates"), Jacobian correction is 0
453 jd = 0.0;
454 }
455
456
apply_force(colvarvalue const & force)457 void colvar::dihedral::apply_force(colvarvalue const &force)
458 {
459 if (!group1->noforce)
460 group1->apply_colvar_force(force.real_value);
461
462 if (!group2->noforce)
463 group2->apply_colvar_force(force.real_value);
464
465 if (!group3->noforce)
466 group3->apply_colvar_force(force.real_value);
467
468 if (!group4->noforce)
469 group4->apply_colvar_force(force.real_value);
470 }
471
472
473 // metrics functions for cvc implementations with a periodicity
474
dist2(colvarvalue const & x1,colvarvalue const & x2) const475 cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
476 colvarvalue const &x2) const
477 {
478 cvm::real diff = x1.real_value - x2.real_value;
479 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
480 return diff * diff;
481 }
482
483
dist2_lgrad(colvarvalue const & x1,colvarvalue const & x2) const484 colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
485 colvarvalue const &x2) const
486 {
487 cvm::real diff = x1.real_value - x2.real_value;
488 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
489 return 2.0 * diff;
490 }
491
492
dist2_rgrad(colvarvalue const & x1,colvarvalue const & x2) const493 colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
494 colvarvalue const &x2) const
495 {
496 cvm::real diff = x1.real_value - x2.real_value;
497 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
498 return (-2.0) * diff;
499 }
500
501
wrap(colvarvalue & x_unwrapped) const502 void colvar::dihedral::wrap(colvarvalue &x_unwrapped) const
503 {
504 if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
505 x_unwrapped.real_value -= 360.0;
506 return;
507 }
508
509 if ((x_unwrapped.real_value - wrap_center) < -180.0) {
510 x_unwrapped.real_value += 360.0;
511 return;
512 }
513 }
514
515
polar_theta(std::string const & conf)516 colvar::polar_theta::polar_theta(std::string const &conf)
517 : cvc(conf)
518 {
519 function_type = "polar_theta";
520 enable(f_cvc_com_based);
521
522 atoms = parse_group(conf, "atoms");
523 init_total_force_params(conf);
524 x.type(colvarvalue::type_scalar);
525 }
526
527
polar_theta()528 colvar::polar_theta::polar_theta()
529 {
530 function_type = "polar_theta";
531 x.type(colvarvalue::type_scalar);
532 }
533
534
calc_value()535 void colvar::polar_theta::calc_value()
536 {
537 cvm::rvector pos = atoms->center_of_mass();
538 r = atoms->center_of_mass().norm();
539 // Internal values of theta and phi are radians
540 theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.;
541 phi = cvm::atan2(pos.y, pos.x);
542 x.real_value = (180.0/PI) * theta;
543 }
544
545
calc_gradients()546 void colvar::polar_theta::calc_gradients()
547 {
548 if (r == 0.)
549 atoms->set_weighted_gradient(cvm::rvector(0., 0., 0.));
550 else
551 atoms->set_weighted_gradient(cvm::rvector(
552 (180.0/PI) * cvm::cos(theta) * cvm::cos(phi) / r,
553 (180.0/PI) * cvm::cos(theta) * cvm::sin(phi) / r,
554 (180.0/PI) * -cvm::sin(theta) / r));
555 }
556
557
apply_force(colvarvalue const & force)558 void colvar::polar_theta::apply_force(colvarvalue const &force)
559 {
560 if (!atoms->noforce)
561 atoms->apply_colvar_force(force.real_value);
562 }
563
564
simple_scalar_dist_functions(polar_theta)565 simple_scalar_dist_functions(polar_theta)
566
567
568 colvar::polar_phi::polar_phi(std::string const &conf)
569 : cvc(conf)
570 {
571 function_type = "polar_phi";
572 period = 360.0;
573 enable(f_cvc_com_based);
574
575 atoms = parse_group(conf, "atoms");
576 init_total_force_params(conf);
577 x.type(colvarvalue::type_scalar);
578 }
579
580
polar_phi()581 colvar::polar_phi::polar_phi()
582 {
583 function_type = "polar_phi";
584 period = 360.0;
585 x.type(colvarvalue::type_scalar);
586 }
587
588
calc_value()589 void colvar::polar_phi::calc_value()
590 {
591 cvm::rvector pos = atoms->center_of_mass();
592 r = atoms->center_of_mass().norm();
593 // Internal values of theta and phi are radians
594 theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.;
595 phi = cvm::atan2(pos.y, pos.x);
596 x.real_value = (180.0/PI) * phi;
597 }
598
599
calc_gradients()600 void colvar::polar_phi::calc_gradients()
601 {
602 atoms->set_weighted_gradient(cvm::rvector(
603 (180.0/PI) * -cvm::sin(phi) / (r*cvm::sin(theta)),
604 (180.0/PI) * cvm::cos(phi) / (r*cvm::sin(theta)),
605 0.));
606 }
607
608
apply_force(colvarvalue const & force)609 void colvar::polar_phi::apply_force(colvarvalue const &force)
610 {
611 if (!atoms->noforce)
612 atoms->apply_colvar_force(force.real_value);
613 }
614
615
616 // Same as dihedral, for polar_phi
617
dist2(colvarvalue const & x1,colvarvalue const & x2) const618 cvm::real colvar::polar_phi::dist2(colvarvalue const &x1,
619 colvarvalue const &x2) const
620 {
621 cvm::real diff = x1.real_value - x2.real_value;
622 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
623 return diff * diff;
624 }
625
626
dist2_lgrad(colvarvalue const & x1,colvarvalue const & x2) const627 colvarvalue colvar::polar_phi::dist2_lgrad(colvarvalue const &x1,
628 colvarvalue const &x2) const
629 {
630 cvm::real diff = x1.real_value - x2.real_value;
631 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
632 return 2.0 * diff;
633 }
634
635
dist2_rgrad(colvarvalue const & x1,colvarvalue const & x2) const636 colvarvalue colvar::polar_phi::dist2_rgrad(colvarvalue const &x1,
637 colvarvalue const &x2) const
638 {
639 cvm::real diff = x1.real_value - x2.real_value;
640 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
641 return (-2.0) * diff;
642 }
643
644
wrap(colvarvalue & x_unwrapped) const645 void colvar::polar_phi::wrap(colvarvalue &x_unwrapped) const
646 {
647 if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
648 x_unwrapped.real_value -= 360.0;
649 return;
650 }
651
652 if ((x_unwrapped.real_value - wrap_center) < -180.0) {
653 x_unwrapped.real_value += 360.0;
654 return;
655 }
656
657 return;
658 }
659