Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

colvarcomp_rotations.C

Go to the documentation of this file.
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 "colvarvalue.h"
00012 #include "colvarparse.h"
00013 #include "colvar.h"
00014 #include "colvarcomp.h"
00015 
00016 
00017 
00018 colvar::orientation::orientation(std::string const &conf)
00019   : cvc()
00020 {
00021   set_function_type("orientation");
00022   disable(f_cvc_explicit_gradient);
00023   x.type(colvarvalue::type_quaternion);
00024   colvar::orientation::init(conf);
00025 }
00026 
00027 
00028 int colvar::orientation::init(std::string const &conf)
00029 {
00030   int error_code = cvc::init(conf);
00031 
00032   atoms = parse_group(conf, "atoms");
00033   ref_pos.reserve(atoms->size());
00034 
00035   if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
00036     cvm::log("Using reference positions from input file.\n");
00037     if (ref_pos.size() != atoms->size()) {
00038       return cvm::error("Error: reference positions do not "
00039                         "match the number of requested atoms.\n", COLVARS_INPUT_ERROR);
00040     }
00041   }
00042 
00043   {
00044     std::string file_name;
00045     if (get_keyval(conf, "refPositionsFile", file_name)) {
00046 
00047       std::string file_col;
00048       double file_col_value=0.0;
00049       if (get_keyval(conf, "refPositionsCol", file_col, std::string(""))) {
00050         // use PDB flags if column is provided
00051         bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0);
00052         if (found && file_col_value==0.0) {
00053           return cvm::error("Error: refPositionsColValue, "
00054                             "if provided, must be non-zero.\n", COLVARS_INPUT_ERROR);
00055         }
00056       }
00057 
00058       ref_pos.resize(atoms->size());
00059       cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
00060                        file_col, file_col_value);
00061     }
00062   }
00063 
00064   if (!ref_pos.size()) {
00065     return cvm::error("Error: must define a set of "
00066                       "reference coordinates.\n", COLVARS_INPUT_ERROR);
00067   }
00068 
00069 
00070   cvm::rvector ref_cog(0.0, 0.0, 0.0);
00071   size_t i;
00072   for (i = 0; i < ref_pos.size(); i++) {
00073     ref_cog += ref_pos[i];
00074   }
00075   ref_cog /= cvm::real(ref_pos.size());
00076   cvm::log("Centering the reference coordinates on the origin by subtracting "
00077            "the center of geometry at "+
00078            cvm::to_str(-1.0 * ref_cog)+"; it is "
00079            "assumed that each atom is the closest "
00080            "periodic image to the center of geometry.\n");
00081   for (i = 0; i < ref_pos.size(); i++) {
00082     ref_pos[i] -= ref_cog;
00083   }
00084 
00085   get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0));
00086 
00087   // initialize rot member data
00088   if (!atoms->noforce) {
00089     rot.request_group2_gradients(atoms->size());
00090   }
00091 
00092   return error_code;
00093 }
00094 
00095 
00096 colvar::orientation::orientation()
00097   : cvc()
00098 {
00099   set_function_type("orientation");
00100   disable(f_cvc_explicit_gradient);
00101   x.type(colvarvalue::type_quaternion);
00102 }
00103 
00104 
00105 void colvar::orientation::calc_value()
00106 {
00107   rot.b_debug_gradients = is_enabled(f_cvc_debug_gradient);
00108   atoms_cog = atoms->center_of_geometry();
00109 
00110   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00111 
00112   if ((rot.q).inner(ref_quat) >= 0.0) {
00113     x.quaternion_value = rot.q;
00114   } else {
00115     x.quaternion_value = -1.0 * rot.q;
00116   }
00117 }
00118 
00119 
00120 void colvar::orientation::calc_gradients()
00121 {
00122   // gradients have already been calculated and stored within the
00123   // member object "rot"; we're not using the "grad" member of each
00124   // atom object, because it only can represent the gradient of a
00125   // scalar colvar
00126 }
00127 
00128 
00129 void colvar::orientation::apply_force(colvarvalue const &force)
00130 {
00131   cvm::quaternion const &FQ = force.quaternion_value;
00132 
00133   if (!atoms->noforce) {
00134     for (size_t ia = 0; ia < atoms->size(); ia++) {
00135       for (size_t i = 0; i < 4; i++) {
00136         (*atoms)[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]);
00137       }
00138     }
00139   }
00140 }
00141 
00142 
00143 cvm::real colvar::orientation::dist2(colvarvalue const &x1,
00144                                      colvarvalue const &x2) const
00145 {
00146   return x1.quaternion_value.dist2(x2);
00147 }
00148 
00149 
00150 colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
00151                                              colvarvalue const &x2) const
00152 {
00153   return x1.quaternion_value.dist2_grad(x2);
00154 }
00155 
00156 
00157 colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
00158                                              colvarvalue const &x2) const
00159 {
00160   return x2.quaternion_value.dist2_grad(x1);
00161 }
00162 
00163 
00164 
00165 colvar::orientation_angle::orientation_angle(std::string const &conf)
00166   : orientation()
00167 {
00168   set_function_type("orientationAngle");
00169   init_as_angle();
00170   enable(f_cvc_explicit_gradient);
00171   orientation_angle::init(conf);
00172 }
00173 
00174 
00175 int colvar::orientation_angle::init(std::string const &conf)
00176 {
00177   return orientation::init(conf);
00178 }
00179 
00180 
00181 void colvar::orientation_angle::calc_value()
00182 {
00183   atoms_cog = atoms->center_of_geometry();
00184 
00185   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00186 
00187   if ((rot.q).q0 >= 0.0) {
00188     x.real_value = (180.0/PI) * 2.0 * cvm::acos((rot.q).q0);
00189   } else {
00190     x.real_value = (180.0/PI) * 2.0 * cvm::acos(-1.0 * (rot.q).q0);
00191   }
00192 }
00193 
00194 
00195 void colvar::orientation_angle::calc_gradients()
00196 {
00197   cvm::real const dxdq0 =
00198     ( ((rot.q).q0 * (rot.q).q0 < 1.0) ?
00199       ((180.0 / PI) * (-2.0) / cvm::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) :
00200       0.0 );
00201 
00202   for (size_t ia = 0; ia < atoms->size(); ia++) {
00203     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
00204   }
00205 }
00206 
00207 
00208 void colvar::orientation_angle::apply_force(colvarvalue const &force)
00209 {
00210   cvm::real const &fw = force.real_value;
00211   if (!atoms->noforce) {
00212     atoms->apply_colvar_force(fw);
00213   }
00214 }
00215 
00216 
00217 simple_scalar_dist_functions(orientation_angle)
00218 
00219 
00220 
00221 colvar::orientation_proj::orientation_proj(std::string const &conf)
00222   : orientation()
00223 {
00224   set_function_type("orientationProj");
00225   enable(f_cvc_explicit_gradient);
00226   x.type(colvarvalue::type_scalar);
00227   init_scalar_boundaries(0.0, 1.0);
00228   orientation_proj::init(conf);
00229 }
00230 
00231 
00232 int colvar::orientation_proj::init(std::string const &conf)
00233 {
00234   return orientation::init(conf);
00235 }
00236 
00237 
00238 void colvar::orientation_proj::calc_value()
00239 {
00240   atoms_cog = atoms->center_of_geometry();
00241   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00242   x.real_value = 2.0 * (rot.q).q0 * (rot.q).q0 - 1.0;
00243 }
00244 
00245 
00246 void colvar::orientation_proj::calc_gradients()
00247 {
00248   cvm::real const dxdq0 = 2.0 * 2.0 * (rot.q).q0;
00249   for (size_t ia = 0; ia < atoms->size(); ia++) {
00250     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
00251   }
00252 }
00253 
00254 
00255 void colvar::orientation_proj::apply_force(colvarvalue const &force)
00256 {
00257   cvm::real const &fw = force.real_value;
00258 
00259   if (!atoms->noforce) {
00260     atoms->apply_colvar_force(fw);
00261   }
00262 }
00263 
00264 
00265 simple_scalar_dist_functions(orientation_proj)
00266 
00267 
00268 
00269 colvar::tilt::tilt(std::string const &conf)
00270   : orientation()
00271 {
00272   set_function_type("tilt");
00273   x.type(colvarvalue::type_scalar);
00274   enable(f_cvc_explicit_gradient);
00275   init_scalar_boundaries(-1.0, 1.0);
00276   tilt::init(conf);
00277 }
00278 
00279 
00280 int colvar::tilt::init(std::string const &conf)
00281 {
00282   int error_code = COLVARS_OK;
00283 
00284   error_code |= orientation::init(conf);
00285 
00286   get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0));
00287   if (axis.norm2() != 1.0) {
00288     axis /= axis.norm();
00289     cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n");
00290   }
00291 
00292   return error_code;
00293 }
00294 
00295 
00296 void colvar::tilt::calc_value()
00297 {
00298   atoms_cog = atoms->center_of_geometry();
00299 
00300   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00301 
00302   x.real_value = rot.cos_theta(axis);
00303 }
00304 
00305 
00306 void colvar::tilt::calc_gradients()
00307 {
00308   cvm::quaternion const dxdq = rot.dcos_theta_dq(axis);
00309 
00310   for (size_t ia = 0; ia < atoms->size(); ia++) {
00311     (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
00312     for (size_t iq = 0; iq < 4; iq++) {
00313       (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
00314     }
00315   }
00316 }
00317 
00318 
00319 void colvar::tilt::apply_force(colvarvalue const &force)
00320 {
00321   cvm::real const &fw = force.real_value;
00322 
00323   if (!atoms->noforce) {
00324     atoms->apply_colvar_force(fw);
00325   }
00326 }
00327 
00328 
00329 simple_scalar_dist_functions(tilt)
00330 
00331 
00332 
00333 colvar::spin_angle::spin_angle(std::string const &conf)
00334   : orientation()
00335 {
00336   set_function_type("spinAngle");
00337   init_as_periodic_angle();
00338   enable(f_cvc_periodic);
00339   enable(f_cvc_explicit_gradient);
00340   spin_angle::init(conf);
00341 }
00342 
00343 
00344 int colvar::spin_angle::init(std::string const &conf)
00345 {
00346   int error_code = COLVARS_OK;
00347 
00348   error_code |= orientation::init(conf);
00349 
00350   get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0));
00351   if (axis.norm2() != 1.0) {
00352     axis /= axis.norm();
00353     cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n");
00354   }
00355 
00356   return error_code;
00357 }
00358 
00359 
00360 colvar::spin_angle::spin_angle()
00361   : orientation()
00362 {
00363   set_function_type("spinAngle");
00364   period = 360.0;
00365   enable(f_cvc_periodic);
00366   enable(f_cvc_explicit_gradient);
00367   x.type(colvarvalue::type_scalar);
00368 }
00369 
00370 
00371 void colvar::spin_angle::calc_value()
00372 {
00373   atoms_cog = atoms->center_of_geometry();
00374 
00375   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00376 
00377   x.real_value = rot.spin_angle(axis);
00378   this->wrap(x);
00379 }
00380 
00381 
00382 void colvar::spin_angle::calc_gradients()
00383 {
00384   cvm::quaternion const dxdq = rot.dspin_angle_dq(axis);
00385 
00386   for (size_t ia = 0; ia < atoms->size(); ia++) {
00387     (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
00388     for (size_t iq = 0; iq < 4; iq++) {
00389       (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
00390     }
00391   }
00392 }
00393 
00394 
00395 void colvar::spin_angle::apply_force(colvarvalue const &force)
00396 {
00397   cvm::real const &fw = force.real_value;
00398 
00399   if (!atoms->noforce) {
00400     atoms->apply_colvar_force(fw);
00401   }
00402 }
00403 
00404 
00405 cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
00406                                     colvarvalue const &x2) const
00407 {
00408   cvm::real diff = x1.real_value - x2.real_value;
00409   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00410   return diff * diff;
00411 }
00412 
00413 
00414 colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
00415                                             colvarvalue const &x2) const
00416 {
00417   cvm::real diff = x1.real_value - x2.real_value;
00418   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00419   return 2.0 * diff;
00420 }
00421 
00422 
00423 colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
00424                                             colvarvalue const &x2) const
00425 {
00426   cvm::real diff = x1.real_value - x2.real_value;
00427   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00428   return (-2.0) * diff;
00429 }
00430 
00431 
00432 void colvar::spin_angle::wrap(colvarvalue &x_unwrapped) const
00433 {
00434   if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
00435     x_unwrapped.real_value -= 360.0;
00436     return;
00437   }
00438 
00439   if ((x_unwrapped.real_value - wrap_center) < -180.0) {
00440     x_unwrapped.real_value += 360.0;
00441     return;
00442   }
00443 
00444   return;
00445 }
00446 
00447 
00448 colvar::euler_phi::euler_phi(std::string const &conf)
00449   : orientation()
00450 {
00451   set_function_type("eulerPhi");
00452   init_as_periodic_angle();
00453   enable(f_cvc_explicit_gradient);
00454   euler_phi::init(conf);
00455 }
00456 
00457 
00458 colvar::euler_phi::euler_phi()
00459   : orientation()
00460 {
00461   set_function_type("eulerPhi");
00462   init_as_periodic_angle();
00463   enable(f_cvc_explicit_gradient);
00464 }
00465 
00466 
00467 int colvar::euler_phi::init(std::string const &conf)
00468 {
00469   int error_code = COLVARS_OK;
00470   error_code |= orientation::init(conf);
00471   return error_code;
00472 }
00473 
00474 
00475 void colvar::euler_phi::calc_value()
00476 {
00477   atoms_cog = atoms->center_of_geometry();
00478 
00479   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00480 
00481   const cvm::real& q0 = rot.q.q0;
00482   const cvm::real& q1 = rot.q.q1;
00483   const cvm::real& q2 = rot.q.q2;
00484   const cvm::real& q3 = rot.q.q3;
00485   const cvm::real tmp_y = 2 * (q0 * q1 + q2 * q3);
00486   const cvm::real tmp_x = 1 - 2 * (q1 * q1 + q2 * q2);
00487   x.real_value = cvm::atan2(tmp_y, tmp_x) * (180.0/PI);
00488 }
00489 
00490 
00491 void colvar::euler_phi::calc_gradients()
00492 {
00493   const cvm::real& q0 = rot.q.q0;
00494   const cvm::real& q1 = rot.q.q1;
00495   const cvm::real& q2 = rot.q.q2;
00496   const cvm::real& q3 = rot.q.q3;
00497   const cvm::real denominator = (2 * q0 * q1 + 2 * q2 * q3) * (2 * q0 * q1 + 2 * q2 * q3) + (-2 * q1 * q1 - 2 * q2 * q2 + 1) * (-2 * q1 * q1 - 2 * q2 * q2 + 1);
00498   const cvm::real dxdq0 = (180.0/PI) * 2 * q1 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) / denominator;
00499   const cvm::real dxdq1 = (180.0/PI) * (2 * q0 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) - 4 * q1 * (-2 * q0 * q1 - 2 * q2 * q3)) / denominator;
00500   const cvm::real dxdq2 = (180.0/PI) * (-4 * q2 * (-2 * q0 * q1 - 2 * q2 * q3) + 2 * q3 * (-2 * q1 * q1 - 2 * q2 * q2 + 1)) / denominator;
00501   const cvm::real dxdq3 = (180.0/PI) * 2 * q2 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) / denominator;
00502   for (size_t ia = 0; ia < atoms->size(); ia++) {
00503     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) +
00504                         (dxdq1 * (rot.dQ0_2[ia])[1]) +
00505                         (dxdq2 * (rot.dQ0_2[ia])[2]) +
00506                         (dxdq3 * (rot.dQ0_2[ia])[3]);
00507   }
00508 }
00509 
00510 
00511 void colvar::euler_phi::apply_force(colvarvalue const &force)
00512 {
00513   cvm::real const &fw = force.real_value;
00514   if (!atoms->noforce) {
00515     atoms->apply_colvar_force(fw);
00516   }
00517 }
00518 
00519 
00520 cvm::real colvar::euler_phi::dist2(colvarvalue const &x1,
00521                                   colvarvalue const &x2) const
00522 {
00523   cvm::real diff = x1.real_value - x2.real_value;
00524   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00525   return diff * diff;
00526 }
00527 
00528 
00529 colvarvalue colvar::euler_phi::dist2_lgrad(colvarvalue const &x1,
00530                                           colvarvalue const &x2) const
00531 {
00532   cvm::real diff = x1.real_value - x2.real_value;
00533   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00534   return 2.0 * diff;
00535 }
00536 
00537 
00538 colvarvalue colvar::euler_phi::dist2_rgrad(colvarvalue const &x1,
00539                                           colvarvalue const &x2) const
00540 {
00541   cvm::real diff = x1.real_value - x2.real_value;
00542   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00543   return (-2.0) * diff;
00544 }
00545 
00546 
00547 void colvar::euler_phi::wrap(colvarvalue &x_unwrapped) const
00548 {
00549   if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
00550     x_unwrapped.real_value -= 360.0;
00551     return;
00552   }
00553 
00554   if ((x_unwrapped.real_value - wrap_center) < -180.0) {
00555     x_unwrapped.real_value += 360.0;
00556     return;
00557   }
00558 
00559   return;
00560 }
00561 
00562 
00563 colvar::euler_psi::euler_psi(std::string const &conf)
00564   : orientation()
00565 {
00566   set_function_type("eulerPsi");
00567   init_as_periodic_angle();
00568   enable(f_cvc_explicit_gradient);
00569   euler_psi::init(conf);
00570 }
00571 
00572 
00573 colvar::euler_psi::euler_psi()
00574   : orientation()
00575 {
00576   set_function_type("eulerPsi");
00577   init_as_periodic_angle();
00578   enable(f_cvc_explicit_gradient);
00579 }
00580 
00581 
00582 int colvar::euler_psi::init(std::string const &conf)
00583 {
00584   int error_code = COLVARS_OK;
00585   error_code |= orientation::init(conf);
00586   return error_code;
00587 }
00588 
00589 
00590 void colvar::euler_psi::calc_value()
00591 {
00592   atoms_cog = atoms->center_of_geometry();
00593 
00594   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00595 
00596   const cvm::real& q0 = rot.q.q0;
00597   const cvm::real& q1 = rot.q.q1;
00598   const cvm::real& q2 = rot.q.q2;
00599   const cvm::real& q3 = rot.q.q3;
00600   const cvm::real tmp_y = 2 * (q0 * q3 + q1 * q2);
00601   const cvm::real tmp_x = 1 - 2 * (q2 * q2 + q3 * q3);
00602   x.real_value = cvm::atan2(tmp_y, tmp_x) * (180.0/PI);
00603 }
00604 
00605 
00606 void colvar::euler_psi::calc_gradients()
00607 {
00608   const cvm::real& q0 = rot.q.q0;
00609   const cvm::real& q1 = rot.q.q1;
00610   const cvm::real& q2 = rot.q.q2;
00611   const cvm::real& q3 = rot.q.q3;
00612   const cvm::real denominator = (2 * q0 * q3 + 2 * q1 * q2) * (2 * q0 * q3 + 2 * q1 * q2) + (-2 * q2 * q2 - 2 * q3 * q3 + 1) * (-2 * q2 * q2 - 2 * q3 * q3 + 1);
00613   const cvm::real dxdq0 = (180.0/PI) * 2 * q3 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) / denominator;
00614   const cvm::real dxdq1 = (180.0/PI) * 2 * q2 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) / denominator;
00615   const cvm::real dxdq2 = (180.0/PI) * (2 * q1 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) - 4 * q2 * (-2 * q0 * q3 - 2 * q1 * q2)) / denominator;
00616   const cvm::real dxdq3 = (180.0/PI) * (2 * q0 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) - 4 * q3 * (-2 * q0 * q3 - 2 * q1 * q2)) / denominator;
00617   for (size_t ia = 0; ia < atoms->size(); ia++) {
00618     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) +
00619                         (dxdq1 * (rot.dQ0_2[ia])[1]) +
00620                         (dxdq2 * (rot.dQ0_2[ia])[2]) +
00621                         (dxdq3 * (rot.dQ0_2[ia])[3]);
00622   }
00623 }
00624 
00625 
00626 void colvar::euler_psi::apply_force(colvarvalue const &force)
00627 {
00628   cvm::real const &fw = force.real_value;
00629   if (!atoms->noforce) {
00630     atoms->apply_colvar_force(fw);
00631   }
00632 }
00633 
00634 
00635 cvm::real colvar::euler_psi::dist2(colvarvalue const &x1,
00636                                   colvarvalue const &x2) const
00637 {
00638   cvm::real diff = x1.real_value - x2.real_value;
00639   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00640   return diff * diff;
00641 }
00642 
00643 
00644 colvarvalue colvar::euler_psi::dist2_lgrad(colvarvalue const &x1,
00645                                           colvarvalue const &x2) const
00646 {
00647   cvm::real diff = x1.real_value - x2.real_value;
00648   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00649   return 2.0 * diff;
00650 }
00651 
00652 
00653 colvarvalue colvar::euler_psi::dist2_rgrad(colvarvalue const &x1,
00654                                           colvarvalue const &x2) const
00655 {
00656   cvm::real diff = x1.real_value - x2.real_value;
00657   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
00658   return (-2.0) * diff;
00659 }
00660 
00661 
00662 void colvar::euler_psi::wrap(colvarvalue &x_unwrapped) const
00663 {
00664   if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
00665     x_unwrapped.real_value -= 360.0;
00666     return;
00667   }
00668 
00669   if ((x_unwrapped.real_value - wrap_center) < -180.0) {
00670     x_unwrapped.real_value += 360.0;
00671     return;
00672   }
00673 
00674   return;
00675 }
00676 
00677 
00678 colvar::euler_theta::euler_theta(std::string const &conf)
00679   : orientation()
00680 {
00681   set_function_type("eulerTheta");
00682   init_as_angle();
00683   enable(f_cvc_explicit_gradient);
00684   euler_theta::init(conf);
00685 }
00686 
00687 
00688 colvar::euler_theta::euler_theta()
00689   : orientation()
00690 {
00691   set_function_type("eulerTheta");
00692   init_as_angle();
00693   enable(f_cvc_explicit_gradient);
00694 }
00695 
00696 
00697 int colvar::euler_theta::init(std::string const &conf)
00698 {
00699   int error_code = COLVARS_OK;
00700   error_code |= orientation::init(conf);
00701   return error_code;
00702 }
00703 
00704 
00705 void colvar::euler_theta::calc_value()
00706 {
00707   atoms_cog = atoms->center_of_geometry();
00708 
00709   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
00710 
00711   const cvm::real& q0 = rot.q.q0;
00712   const cvm::real& q1 = rot.q.q1;
00713   const cvm::real& q2 = rot.q.q2;
00714   const cvm::real& q3 = rot.q.q3;
00715   x.real_value = cvm::asin(2 * (q0 * q2 - q3 * q1)) * (180.0/PI);
00716 }
00717 
00718 
00719 void colvar::euler_theta::calc_gradients()
00720 {
00721   const cvm::real& q0 = rot.q.q0;
00722   const cvm::real& q1 = rot.q.q1;
00723   const cvm::real& q2 = rot.q.q2;
00724   const cvm::real& q3 = rot.q.q3;
00725   const cvm::real denominator = cvm::sqrt(1 - (2 * q0 * q2 - 2 * q1 * q3) * (2 * q0 * q2 - 2 * q1 * q3));
00726   const cvm::real dxdq0 = (180.0/PI) * 2 * q2 / denominator;
00727   const cvm::real dxdq1 = (180.0/PI) * -2 * q3 / denominator;
00728   const cvm::real dxdq2 = (180.0/PI) * 2 * q0 / denominator;
00729   const cvm::real dxdq3 = (180.0/PI) * -2 * q1 / denominator;
00730   for (size_t ia = 0; ia < atoms->size(); ia++) {
00731     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) +
00732                         (dxdq1 * (rot.dQ0_2[ia])[1]) +
00733                         (dxdq2 * (rot.dQ0_2[ia])[2]) +
00734                         (dxdq3 * (rot.dQ0_2[ia])[3]);
00735   }
00736 }
00737 
00738 
00739 void colvar::euler_theta::apply_force(colvarvalue const &force)
00740 {
00741   cvm::real const &fw = force.real_value;
00742   if (!atoms->noforce) {
00743     atoms->apply_colvar_force(fw);
00744   }
00745 }
00746 
00747 
00748 cvm::real colvar::euler_theta::dist2(colvarvalue const &x1,
00749                                      colvarvalue const &x2) const
00750 {
00751   // theta angle is not periodic
00752   return cvc::dist2(x1, x2);
00753 }
00754 
00755 
00756 colvarvalue colvar::euler_theta::dist2_lgrad(colvarvalue const &x1,
00757                                              colvarvalue const &x2) const
00758 {
00759   // theta angle is not periodic
00760   return cvc::dist2_lgrad(x1, x2);
00761 }
00762 
00763 
00764 colvarvalue colvar::euler_theta::dist2_rgrad(colvarvalue const &x1,
00765                                              colvarvalue const &x2) const
00766 {
00767   // theta angle is not periodic
00768   return cvc::dist2_rgrad(x1, x2);
00769 }

Generated on Fri Oct 4 02:43:42 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002