00001 // -*- c++ -*- 00002 00003 // This file is part of the Collective Variables module (Colvars). 00004 // The original version of Colvars and its updates are located at: 00005 // https://github.com/Colvars/colvars 00006 // Please update all Colvars source files before making any changes. 00007 // If you wish to distribute your changes, please submit them to the 00008 // Colvars repository at GitHub. 00009 00010 #include "colvarmodule.h" 00011 #include "colvar.h" 00012 #include "colvarcomp.h" 00013 00014 00015 colvar::angle::angle(std::string const &conf) 00016 : cvc(conf) 00017 { 00018 set_function_type("angle"); 00019 init_as_angle(); 00020 00021 provide(f_cvc_inv_gradient); 00022 provide(f_cvc_Jacobian); 00023 enable(f_cvc_com_based); 00024 00025 group1 = parse_group(conf, "group1"); 00026 group2 = parse_group(conf, "group2"); 00027 group3 = parse_group(conf, "group3"); 00028 00029 init_total_force_params(conf); 00030 } 00031 00032 00033 colvar::angle::angle(cvm::atom const &a1, 00034 cvm::atom const &a2, 00035 cvm::atom const &a3) 00036 { 00037 set_function_type("angle"); 00038 init_as_angle(); 00039 00040 provide(f_cvc_inv_gradient); 00041 provide(f_cvc_Jacobian); 00042 enable(f_cvc_com_based); 00043 00044 group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1)); 00045 group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2)); 00046 group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3)); 00047 register_atom_group(group1); 00048 register_atom_group(group2); 00049 register_atom_group(group3); 00050 } 00051 00052 00053 void colvar::angle::calc_value() 00054 { 00055 cvm::atom_pos const g1_pos = group1->center_of_mass(); 00056 cvm::atom_pos const g2_pos = group2->center_of_mass(); 00057 cvm::atom_pos const g3_pos = group3->center_of_mass(); 00058 00059 r21 = is_enabled(f_cvc_pbc_minimum_image) ? 00060 cvm::position_distance(g2_pos, g1_pos) : 00061 g1_pos - g2_pos; 00062 r21l = r21.norm(); 00063 r23 = is_enabled(f_cvc_pbc_minimum_image) ? 00064 cvm::position_distance(g2_pos, g3_pos) : 00065 g3_pos - g2_pos; 00066 r23l = r23.norm(); 00067 00068 cvm::real const cos_theta = (r21*r23)/(r21l*r23l); 00069 00070 x.real_value = (180.0/PI) * cvm::acos(cos_theta); 00071 } 00072 00073 00074 void colvar::angle::calc_gradients() 00075 { 00076 cvm::real const cos_theta = (r21*r23)/(r21l*r23l); 00077 cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta); 00078 00079 dxdr1 = (180.0/PI) * dxdcos * 00080 (1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l ); 00081 00082 dxdr3 = (180.0/PI) * dxdcos * 00083 (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l ); 00084 00085 group1->set_weighted_gradient(dxdr1); 00086 group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0)); 00087 group3->set_weighted_gradient(dxdr3); 00088 } 00089 00090 00091 void colvar::angle::calc_force_invgrads() 00092 { 00093 // This uses a force measurement on groups 1 and 3 only 00094 // to keep in line with the implicit variable change used to 00095 // evaluate the Jacobian term (essentially polar coordinates 00096 // centered on group2, which means group2 is kept fixed 00097 // when propagating changes in the angle) 00098 00099 if (is_enabled(f_cvc_one_site_total_force)) { 00100 group1->read_total_forces(); 00101 cvm::real norm_fact = 1.0 / dxdr1.norm2(); 00102 ft.real_value = norm_fact * dxdr1 * group1->total_force(); 00103 } else { 00104 group1->read_total_forces(); 00105 group3->read_total_forces(); 00106 cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2()); 00107 ft.real_value = norm_fact * ( dxdr1 * group1->total_force() 00108 + dxdr3 * group3->total_force()); 00109 } 00110 return; 00111 } 00112 00113 00114 void colvar::angle::calc_Jacobian_derivative() 00115 { 00116 // det(J) = (2 pi) r^2 * sin(theta) 00117 // hence Jd = cot(theta) 00118 const cvm::real theta = x.real_value * PI / 180.0; 00119 jd = PI / 180.0 * (theta != 0.0 ? cvm::cos(theta) / cvm::sin(theta) : 0.0); 00120 } 00121 00122 00123 void colvar::angle::apply_force(colvarvalue const &force) 00124 { 00125 if (!group1->noforce) 00126 group1->apply_colvar_force(force.real_value); 00127 00128 if (!group2->noforce) 00129 group2->apply_colvar_force(force.real_value); 00130 00131 if (!group3->noforce) 00132 group3->apply_colvar_force(force.real_value); 00133 } 00134 00135 00136 simple_scalar_dist_functions(angle) 00137 00138 00139 00140 colvar::dipole_angle::dipole_angle(std::string const &conf) 00141 : cvc(conf) 00142 { 00143 set_function_type("dipoleAngle"); 00144 init_as_angle(); 00145 00146 group1 = parse_group(conf, "group1"); 00147 group2 = parse_group(conf, "group2"); 00148 group3 = parse_group(conf, "group3"); 00149 00150 init_total_force_params(conf); 00151 } 00152 00153 00154 colvar::dipole_angle::dipole_angle(cvm::atom const &a1, 00155 cvm::atom const &a2, 00156 cvm::atom const &a3) 00157 { 00158 set_function_type("dipoleAngle"); 00159 init_as_angle(); 00160 00161 group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1)); 00162 group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2)); 00163 group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3)); 00164 register_atom_group(group1); 00165 register_atom_group(group2); 00166 register_atom_group(group3); 00167 } 00168 00169 00170 colvar::dipole_angle::dipole_angle() 00171 { 00172 set_function_type("dipoleAngle"); 00173 init_as_angle(); 00174 } 00175 00176 00177 void colvar::dipole_angle::calc_value() 00178 { 00179 cvm::atom_pos const g1_pos = group1->center_of_mass(); 00180 cvm::atom_pos const g2_pos = group2->center_of_mass(); 00181 cvm::atom_pos const g3_pos = group3->center_of_mass(); 00182 00183 group1->calc_dipole(g1_pos); 00184 00185 r21 = group1->dipole(); 00186 r21l = r21.norm(); 00187 r23 = is_enabled(f_cvc_pbc_minimum_image) ? 00188 cvm::position_distance(g2_pos, g3_pos) : 00189 g3_pos - g2_pos; 00190 r23l = r23.norm(); 00191 00192 cvm::real const cos_theta = (r21*r23)/(r21l*r23l); 00193 00194 x.real_value = (180.0/PI) * cvm::acos(cos_theta); 00195 } 00196 00197 //to be implemented 00198 //void colvar::dipole_angle::calc_force_invgrads(){} 00199 //void colvar::dipole_angle::calc_Jacobian_derivative(){} 00200 00201 void colvar::dipole_angle::calc_gradients() 00202 { 00203 cvm::real const cos_theta = (r21*r23)/(r21l*r23l); 00204 cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta); 00205 00206 dxdr1 = (180.0/PI) * dxdcos * 00207 (1.0/r21l)* (r23/r23l + (-1.0) * cos_theta * r21/r21l ); 00208 00209 dxdr3 = (180.0/PI) * dxdcos * 00210 (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l ); 00211 00212 //this auxiliar variables are to avoid numerical errors inside "for" 00213 double aux1 = group1->total_charge/group1->total_mass; 00214 // double aux2 = group2->total_charge/group2->total_mass; 00215 // double aux3 = group3->total_charge/group3->total_mass; 00216 00217 size_t i; 00218 for (i = 0; i < group1->size(); i++) { 00219 (*group1)[i].grad =((*group1)[i].charge + (-1)* (*group1)[i].mass * aux1) * (dxdr1); 00220 } 00221 00222 for (i = 0; i < group2->size(); i++) { 00223 (*group2)[i].grad = ((*group2)[i].mass/group2->total_mass)* dxdr3 * (-1.0); 00224 } 00225 00226 for (i = 0; i < group3->size(); i++) { 00227 (*group3)[i].grad =((*group3)[i].mass/group3->total_mass) * (dxdr3); 00228 } 00229 } 00230 00231 00232 void colvar::dipole_angle::apply_force(colvarvalue const &force) 00233 { 00234 if (!group1->noforce) 00235 group1->apply_colvar_force(force.real_value); 00236 00237 if (!group2->noforce) 00238 group2->apply_colvar_force(force.real_value); 00239 00240 if (!group3->noforce) 00241 group3->apply_colvar_force(force.real_value); 00242 } 00243 00244 00245 simple_scalar_dist_functions(dipole_angle) 00246 00247 00248 00249 colvar::dihedral::dihedral(std::string const &conf) 00250 : cvc(conf) 00251 { 00252 set_function_type("dihedral"); 00253 init_as_periodic_angle(); 00254 provide(f_cvc_inv_gradient); 00255 provide(f_cvc_Jacobian); 00256 enable(f_cvc_com_based); 00257 00258 group1 = parse_group(conf, "group1"); 00259 group2 = parse_group(conf, "group2"); 00260 group3 = parse_group(conf, "group3"); 00261 group4 = parse_group(conf, "group4"); 00262 00263 init_total_force_params(conf); 00264 } 00265 00266 00267 colvar::dihedral::dihedral(cvm::atom const &a1, 00268 cvm::atom const &a2, 00269 cvm::atom const &a3, 00270 cvm::atom const &a4) 00271 { 00272 set_function_type("dihedral"); 00273 init_as_periodic_angle(); 00274 provide(f_cvc_inv_gradient); 00275 provide(f_cvc_Jacobian); 00276 enable(f_cvc_com_based); 00277 00278 b_1site_force = false; 00279 00280 group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1)); 00281 group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2)); 00282 group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3)); 00283 group4 = new cvm::atom_group(std::vector<cvm::atom>(1, a4)); 00284 register_atom_group(group1); 00285 register_atom_group(group2); 00286 register_atom_group(group3); 00287 register_atom_group(group4); 00288 } 00289 00290 00291 colvar::dihedral::dihedral() 00292 { 00293 set_function_type("dihedral"); 00294 init_as_periodic_angle(); 00295 enable(f_cvc_periodic); 00296 provide(f_cvc_inv_gradient); 00297 provide(f_cvc_Jacobian); 00298 } 00299 00300 00301 void colvar::dihedral::calc_value() 00302 { 00303 cvm::atom_pos const g1_pos = group1->center_of_mass(); 00304 cvm::atom_pos const g2_pos = group2->center_of_mass(); 00305 cvm::atom_pos const g3_pos = group3->center_of_mass(); 00306 cvm::atom_pos const g4_pos = group4->center_of_mass(); 00307 00308 // Usual sign convention: r12 = r2 - r1 00309 r12 = is_enabled(f_cvc_pbc_minimum_image) ? 00310 cvm::position_distance(g1_pos, g2_pos) : 00311 g2_pos - g1_pos; 00312 r23 = is_enabled(f_cvc_pbc_minimum_image) ? 00313 cvm::position_distance(g2_pos, g3_pos) : 00314 g3_pos - g2_pos; 00315 r34 = is_enabled(f_cvc_pbc_minimum_image) ? 00316 cvm::position_distance(g3_pos, g4_pos) : 00317 g4_pos - g3_pos; 00318 00319 cvm::rvector const n1 = cvm::rvector::outer(r12, r23); 00320 cvm::rvector const n2 = cvm::rvector::outer(r23, r34); 00321 00322 cvm::real const cos_phi = n1 * n2; 00323 cvm::real const sin_phi = n1 * r34 * r23.norm(); 00324 00325 x.real_value = (180.0/PI) * cvm::atan2(sin_phi, cos_phi); 00326 this->wrap(x); 00327 } 00328 00329 00330 void colvar::dihedral::calc_gradients() 00331 { 00332 cvm::rvector A = cvm::rvector::outer(r12, r23); 00333 cvm::real rA = A.norm(); 00334 cvm::rvector B = cvm::rvector::outer(r23, r34); 00335 cvm::real rB = B.norm(); 00336 cvm::rvector C = cvm::rvector::outer(r23, A); 00337 cvm::real rC = C.norm(); 00338 00339 cvm::real const cos_phi = (A*B)/(rA*rB); 00340 cvm::real const sin_phi = (C*B)/(rC*rB); 00341 00342 cvm::rvector f1, f2, f3; 00343 00344 rB = 1.0/rB; 00345 B *= rB; 00346 00347 if (cvm::fabs(sin_phi) > 0.1) { 00348 rA = 1.0/rA; 00349 A *= rA; 00350 cvm::rvector const dcosdA = rA*(cos_phi*A-B); 00351 cvm::rvector const dcosdB = rB*(cos_phi*B-A); 00352 rA = 1.0; 00353 00354 cvm::real const K = (1.0/sin_phi) * (180.0/PI); 00355 00356 f1 = K * cvm::rvector::outer(r23, dcosdA); 00357 f3 = K * cvm::rvector::outer(dcosdB, r23); 00358 f2 = K * (cvm::rvector::outer(dcosdA, r12) 00359 + cvm::rvector::outer(r34, dcosdB)); 00360 } 00361 else { 00362 rC = 1.0/rC; 00363 C *= rC; 00364 cvm::rvector const dsindC = rC*(sin_phi*C-B); 00365 cvm::rvector const dsindB = rB*(sin_phi*B-C); 00366 rC = 1.0; 00367 00368 cvm::real const K = (-1.0/cos_phi) * (180.0/PI); 00369 00370 f1.x = K*((r23.y*r23.y + r23.z*r23.z)*dsindC.x 00371 - r23.x*r23.y*dsindC.y 00372 - r23.x*r23.z*dsindC.z); 00373 f1.y = K*((r23.z*r23.z + r23.x*r23.x)*dsindC.y 00374 - r23.y*r23.z*dsindC.z 00375 - r23.y*r23.x*dsindC.x); 00376 f1.z = K*((r23.x*r23.x + r23.y*r23.y)*dsindC.z 00377 - r23.z*r23.x*dsindC.x 00378 - r23.z*r23.y*dsindC.y); 00379 00380 f3 = cvm::rvector::outer(dsindB, r23); 00381 f3 *= K; 00382 00383 f2.x = K*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x 00384 +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y 00385 +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z 00386 +dsindB.z*r34.y - dsindB.y*r34.z); 00387 f2.y = K*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y 00388 +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z 00389 +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x 00390 +dsindB.x*r34.z - dsindB.z*r34.x); 00391 f2.z = K*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z 00392 +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x 00393 +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y 00394 +dsindB.y*r34.x - dsindB.x*r34.y); 00395 } 00396 00397 group1->set_weighted_gradient(-f1); 00398 group2->set_weighted_gradient(-f2 + f1); 00399 group3->set_weighted_gradient(-f3 + f2); 00400 group4->set_weighted_gradient(f3); 00401 } 00402 00403 00404 void colvar::dihedral::calc_force_invgrads() 00405 { 00406 cvm::rvector const u12 = r12.unit(); 00407 cvm::rvector const u23 = r23.unit(); 00408 cvm::rvector const u34 = r34.unit(); 00409 00410 cvm::real const d12 = r12.norm(); 00411 cvm::real const d34 = r34.norm(); 00412 00413 cvm::rvector const cross1 = (cvm::rvector::outer(u23, u12)).unit(); 00414 cvm::rvector const cross4 = (cvm::rvector::outer(u23, u34)).unit(); 00415 00416 cvm::real const dot1 = u23 * u12; 00417 cvm::real const dot4 = u23 * u34; 00418 00419 cvm::real const fact1 = d12 * cvm::sqrt(1.0 - dot1 * dot1); 00420 cvm::real const fact4 = d34 * cvm::sqrt(1.0 - dot4 * dot4); 00421 00422 group1->read_total_forces(); 00423 if (is_enabled(f_cvc_one_site_total_force)) { 00424 // This is only measuring the force on group 1 00425 ft.real_value = PI/180.0 * fact1 * (cross1 * group1->total_force()); 00426 } else { 00427 // Default case: use groups 1 and 4 00428 group4->read_total_forces(); 00429 ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1->total_force()) 00430 + fact4 * (cross4 * group4->total_force())); 00431 } 00432 } 00433 00434 00435 void colvar::dihedral::calc_Jacobian_derivative() 00436 { 00437 // With this choice of inverse gradient ("internal coordinates"), Jacobian correction is 0 00438 jd = 0.0; 00439 } 00440 00441 00442 void colvar::dihedral::apply_force(colvarvalue const &force) 00443 { 00444 if (!group1->noforce) 00445 group1->apply_colvar_force(force.real_value); 00446 00447 if (!group2->noforce) 00448 group2->apply_colvar_force(force.real_value); 00449 00450 if (!group3->noforce) 00451 group3->apply_colvar_force(force.real_value); 00452 00453 if (!group4->noforce) 00454 group4->apply_colvar_force(force.real_value); 00455 } 00456 00457 00458 // metrics functions for cvc implementations with a periodicity 00459 00460 cvm::real colvar::dihedral::dist2(colvarvalue const &x1, 00461 colvarvalue const &x2) const 00462 { 00463 cvm::real diff = x1.real_value - x2.real_value; 00464 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); 00465 return diff * diff; 00466 } 00467 00468 00469 colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1, 00470 colvarvalue const &x2) const 00471 { 00472 cvm::real diff = x1.real_value - x2.real_value; 00473 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); 00474 return 2.0 * diff; 00475 } 00476 00477 00478 colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1, 00479 colvarvalue const &x2) const 00480 { 00481 cvm::real diff = x1.real_value - x2.real_value; 00482 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); 00483 return (-2.0) * diff; 00484 } 00485 00486 00487 void colvar::dihedral::wrap(colvarvalue &x_unwrapped) const 00488 { 00489 if ((x_unwrapped.real_value - wrap_center) >= 180.0) { 00490 x_unwrapped.real_value -= 360.0; 00491 return; 00492 } 00493 00494 if ((x_unwrapped.real_value - wrap_center) < -180.0) { 00495 x_unwrapped.real_value += 360.0; 00496 return; 00497 } 00498 } 00499 00500 00501 colvar::polar_theta::polar_theta(std::string const &conf) 00502 : cvc(conf) 00503 { 00504 set_function_type("polarTheta"); 00505 enable(f_cvc_com_based); 00506 00507 atoms = parse_group(conf, "atoms"); 00508 init_total_force_params(conf); 00509 x.type(colvarvalue::type_scalar); 00510 } 00511 00512 00513 colvar::polar_theta::polar_theta() 00514 { 00515 set_function_type("polarTheta"); 00516 x.type(colvarvalue::type_scalar); 00517 } 00518 00519 00520 void colvar::polar_theta::calc_value() 00521 { 00522 cvm::rvector pos = atoms->center_of_mass(); 00523 r = atoms->center_of_mass().norm(); 00524 // Internal values of theta and phi are radians 00525 theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.; 00526 phi = cvm::atan2(pos.y, pos.x); 00527 x.real_value = (180.0/PI) * theta; 00528 } 00529 00530 00531 void colvar::polar_theta::calc_gradients() 00532 { 00533 if (r == 0.) 00534 atoms->set_weighted_gradient(cvm::rvector(0., 0., 0.)); 00535 else 00536 atoms->set_weighted_gradient(cvm::rvector( 00537 (180.0/PI) * cvm::cos(theta) * cvm::cos(phi) / r, 00538 (180.0/PI) * cvm::cos(theta) * cvm::sin(phi) / r, 00539 (180.0/PI) * -cvm::sin(theta) / r)); 00540 } 00541 00542 00543 void colvar::polar_theta::apply_force(colvarvalue const &force) 00544 { 00545 if (!atoms->noforce) 00546 atoms->apply_colvar_force(force.real_value); 00547 } 00548 00549 00550 simple_scalar_dist_functions(polar_theta) 00551 00552 00553 colvar::polar_phi::polar_phi(std::string const &conf) 00554 : cvc(conf) 00555 { 00556 set_function_type("polarPhi"); 00557 init_as_periodic_angle(); 00558 enable(f_cvc_com_based); 00559 00560 atoms = parse_group(conf, "atoms"); 00561 init_total_force_params(conf); 00562 } 00563 00564 00565 colvar::polar_phi::polar_phi() 00566 { 00567 set_function_type("polarPhi"); 00568 init_as_periodic_angle(); 00569 } 00570 00571 00572 void colvar::polar_phi::calc_value() 00573 { 00574 cvm::rvector pos = atoms->center_of_mass(); 00575 r = atoms->center_of_mass().norm(); 00576 // Internal values of theta and phi are radians 00577 theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.; 00578 phi = cvm::atan2(pos.y, pos.x); 00579 x.real_value = (180.0/PI) * phi; 00580 } 00581 00582 00583 void colvar::polar_phi::calc_gradients() 00584 { 00585 atoms->set_weighted_gradient(cvm::rvector( 00586 (180.0/PI) * -cvm::sin(phi) / (r*cvm::sin(theta)), 00587 (180.0/PI) * cvm::cos(phi) / (r*cvm::sin(theta)), 00588 0.)); 00589 } 00590 00591 00592 void colvar::polar_phi::apply_force(colvarvalue const &force) 00593 { 00594 if (!atoms->noforce) 00595 atoms->apply_colvar_force(force.real_value); 00596 } 00597 00598 00599 // Same as dihedral, for polar_phi 00600 00601 cvm::real colvar::polar_phi::dist2(colvarvalue const &x1, 00602 colvarvalue const &x2) const 00603 { 00604 cvm::real diff = x1.real_value - x2.real_value; 00605 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); 00606 return diff * diff; 00607 } 00608 00609 00610 colvarvalue colvar::polar_phi::dist2_lgrad(colvarvalue const &x1, 00611 colvarvalue const &x2) const 00612 { 00613 cvm::real diff = x1.real_value - x2.real_value; 00614 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); 00615 return 2.0 * diff; 00616 } 00617 00618 00619 colvarvalue colvar::polar_phi::dist2_rgrad(colvarvalue const &x1, 00620 colvarvalue const &x2) const 00621 { 00622 cvm::real diff = x1.real_value - x2.real_value; 00623 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); 00624 return (-2.0) * diff; 00625 } 00626 00627 00628 void colvar::polar_phi::wrap(colvarvalue &x_unwrapped) const 00629 { 00630 if ((x_unwrapped.real_value - wrap_center) >= 180.0) { 00631 x_unwrapped.real_value -= 360.0; 00632 return; 00633 } 00634 00635 if ((x_unwrapped.real_value - wrap_center) < -180.0) { 00636 x_unwrapped.real_value += 360.0; 00637 return; 00638 } 00639 00640 return; 00641 }