00001 #include <cmath>
00002
00003 #include "colvarmodule.h"
00004 #include "colvarvalue.h"
00005 #include "colvarparse.h"
00006 #include "colvar.h"
00007 #include "colvarcomp.h"
00008
00009
00010
00011
00012
00013
00014 colvar::distance::distance (std::string const &conf, bool twogroups)
00015 : cvc (conf)
00016 {
00017 function_type = "distance";
00018 b_inverse_gradients = true;
00019 b_Jacobian_derivative = true;
00020 if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) {
00021 cvm::log ("Computing distance using absolute positions (not minimal-image)");
00022 }
00023 if (twogroups && get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
00024 cvm::log ("Computing system force on group 1 only");
00025 }
00026 parse_group (conf, "group1", group1);
00027 atom_groups.push_back (&group1);
00028 if (twogroups) {
00029 parse_group (conf, "group2", group2);
00030 atom_groups.push_back (&group2);
00031 }
00032 x.type (colvarvalue::type_scalar);
00033 }
00034
00035
00036 colvar::distance::distance()
00037 : cvc ()
00038 {
00039 function_type = "distance";
00040 b_inverse_gradients = true;
00041 b_Jacobian_derivative = true;
00042 b_1site_force = false;
00043 x.type (colvarvalue::type_scalar);
00044 }
00045
00046 void colvar::distance::calc_value()
00047 {
00048 if (b_no_PBC) {
00049 dist_v = group2.center_of_mass() - group1.center_of_mass();
00050 } else {
00051 dist_v = cvm::position_distance (group1.center_of_mass(),
00052 group2.center_of_mass());
00053 }
00054 x.real_value = dist_v.norm();
00055 }
00056
00057 void colvar::distance::calc_gradients()
00058 {
00059 cvm::rvector const u = dist_v.unit();
00060 group1.set_weighted_gradient (-1.0 * u);
00061 group2.set_weighted_gradient ( u);
00062 }
00063
00064 void colvar::distance::calc_force_invgrads()
00065 {
00066 group1.read_system_forces();
00067 if ( b_1site_force ) {
00068 ft.real_value = -1.0 * (group1.system_force() * dist_v.unit());
00069 } else {
00070 group2.read_system_forces();
00071 ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit());
00072 }
00073 }
00074
00075 void colvar::distance::calc_Jacobian_derivative()
00076 {
00077 jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;
00078 }
00079
00080 void colvar::distance::apply_force (colvarvalue const &force)
00081 {
00082 if (!group1.noforce)
00083 group1.apply_colvar_force (force.real_value);
00084
00085 if (!group2.noforce)
00086 group2.apply_colvar_force (force.real_value);
00087 }
00088
00089
00090
00091 colvar::distance_vec::distance_vec (std::string const &conf)
00092 : distance (conf)
00093 {
00094 function_type = "distance_vec";
00095 x.type (colvarvalue::type_vector);
00096 }
00097
00098 colvar::distance_vec::distance_vec()
00099 : distance()
00100 {
00101 function_type = "distance_vec";
00102 x.type (colvarvalue::type_vector);
00103 }
00104
00105 void colvar::distance_vec::calc_value()
00106 {
00107 if (b_no_PBC) {
00108 x.rvector_value = group2.center_of_mass() - group1.center_of_mass();
00109 } else {
00110 x.rvector_value = cvm::position_distance (group1.center_of_mass(),
00111 group2.center_of_mass());
00112 }
00113 }
00114
00115 void colvar::distance_vec::calc_gradients()
00116 {
00117
00118
00119 }
00120
00121 void colvar::distance_vec::apply_force (colvarvalue const &force)
00122 {
00123 if (!group1.noforce)
00124 group1.apply_force (-1.0 * force.rvector_value);
00125
00126 if (!group2.noforce)
00127 group2.apply_force ( force.rvector_value);
00128 }
00129
00130
00131
00132 colvar::distance_z::distance_z (std::string const &conf)
00133 : cvc (conf)
00134 {
00135 function_type = "distance_z";
00136 b_inverse_gradients = true;
00137 b_Jacobian_derivative = true;
00138 x.type (colvarvalue::type_scalar);
00139
00140
00141
00142 if (period != 0.0)
00143 b_periodic = true;
00144
00145 if ((wrap_center != 0.0) && (period == 0.0)) {
00146 cvm::fatal_error ("Error: wrapAround was defined in a distanceZ component,"
00147 " but its period has not been set.\n");
00148 }
00149
00150 parse_group (conf, "main", main);
00151 parse_group (conf, "ref", ref1);
00152 atom_groups.push_back (&main);
00153 atom_groups.push_back (&ref1);
00154
00155 parse_group (conf, "ref2", ref2, true);
00156
00157 if (ref2.size()) {
00158 atom_groups.push_back (&ref2);
00159 cvm::log ("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
00160 fixed_axis = false;
00161 if (key_lookup (conf, "axis"))
00162 cvm::log ("Warning: explicit axis definition will be ignored!");
00163 } else {
00164 if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
00165 if (axis.norm2() == 0.0)
00166 cvm::fatal_error ("Axis vector is zero!");
00167 if (axis.norm2() != 1.0) {
00168 axis = axis.unit();
00169 cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
00170 }
00171 }
00172 fixed_axis = true;
00173 }
00174
00175 if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) {
00176 cvm::log ("Computing distance using absolute positions (not minimal-image)");
00177 }
00178 if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
00179 cvm::log ("Computing system force on group \"main\" only");
00180 }
00181 }
00182
00183 colvar::distance_z::distance_z()
00184 {
00185 function_type = "distance_z";
00186 b_inverse_gradients = true;
00187 b_Jacobian_derivative = true;
00188 x.type (colvarvalue::type_scalar);
00189 }
00190
00191 void colvar::distance_z::calc_value()
00192 {
00193 if (fixed_axis) {
00194 if (b_no_PBC) {
00195 dist_v = main.center_of_mass() - ref1.center_of_mass();
00196 } else {
00197 dist_v = cvm::position_distance (ref1.center_of_mass(),
00198 main.center_of_mass());
00199 }
00200 } else {
00201
00202 if (b_no_PBC) {
00203 dist_v = main.center_of_mass() -
00204 (0.5 * (ref1.center_of_mass() + ref2.center_of_mass()));
00205 axis = ref2.center_of_mass() - ref1.center_of_mass();
00206 } else {
00207 dist_v = cvm::position_distance (0.5 * (ref1.center_of_mass() +
00208 ref2.center_of_mass()), main.center_of_mass());
00209 axis = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
00210 }
00211 axis_norm = axis.norm();
00212 axis = axis.unit();
00213 }
00214 x.real_value = axis * dist_v;
00215 this->wrap (x);
00216 }
00217
00218 void colvar::distance_z::calc_gradients()
00219 {
00220 main.set_weighted_gradient ( axis );
00221
00222 if (fixed_axis) {
00223 ref1.set_weighted_gradient (-1.0 * axis);
00224 } else {
00225 if (b_no_PBC) {
00226 ref1.set_weighted_gradient ( 1.0 / axis_norm * (main.center_of_mass() - ref2.center_of_mass() -
00227 x.real_value * axis ));
00228 ref2.set_weighted_gradient ( 1.0 / axis_norm * (ref1.center_of_mass() - main.center_of_mass() +
00229 x.real_value * axis ));
00230 } else {
00231 ref1.set_weighted_gradient ( 1.0 / axis_norm * (
00232 cvm::position_distance (ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis ));
00233 ref2.set_weighted_gradient ( 1.0 / axis_norm * (
00234 cvm::position_distance (main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis ));
00235 }
00236 }
00237 }
00238
00239 void colvar::distance_z::calc_force_invgrads()
00240 {
00241 main.read_system_forces();
00242
00243 if (fixed_axis && !b_1site_force) {
00244 ref1.read_system_forces();
00245 ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis);
00246 } else {
00247 ft.real_value = main.system_force() * axis;
00248 }
00249 }
00250
00251 void colvar::distance_z::calc_Jacobian_derivative()
00252 {
00253 jd.real_value = 0.0;
00254 }
00255
00256 void colvar::distance_z::apply_force (colvarvalue const &force)
00257 {
00258 if (!ref1.noforce)
00259 ref1.apply_colvar_force (force.real_value);
00260
00261 if (ref2.size() && !ref2.noforce)
00262 ref2.apply_colvar_force (force.real_value);
00263
00264 if (!main.noforce)
00265 main.apply_colvar_force (force.real_value);
00266 }
00267
00268
00269
00270 colvar::distance_xy::distance_xy (std::string const &conf)
00271 : distance_z (conf)
00272 {
00273 function_type = "distance_xy";
00274 b_inverse_gradients = true;
00275 b_Jacobian_derivative = true;
00276 x.type (colvarvalue::type_scalar);
00277 }
00278
00279 colvar::distance_xy::distance_xy()
00280 : distance_z()
00281 {
00282 function_type = "distance_xy";
00283 b_inverse_gradients = true;
00284 b_Jacobian_derivative = true;
00285 x.type (colvarvalue::type_scalar);
00286 }
00287
00288 void colvar::distance_xy::calc_value()
00289 {
00290 if (b_no_PBC) {
00291 dist_v = main.center_of_mass() - ref1.center_of_mass();
00292 } else {
00293 dist_v = cvm::position_distance (ref1.center_of_mass(),
00294 main.center_of_mass());
00295 }
00296 if (!fixed_axis) {
00297 if (b_no_PBC) {
00298 v12 = ref2.center_of_mass() - ref1.center_of_mass();
00299 } else {
00300 v12 = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
00301 }
00302 axis_norm = v12.norm();
00303 axis = v12.unit();
00304 }
00305
00306 dist_v_ortho = dist_v - (dist_v * axis) * axis;
00307 x.real_value = dist_v_ortho.norm();
00308 }
00309
00310 void colvar::distance_xy::calc_gradients()
00311 {
00312
00313
00314 cvm::real A;
00315 cvm::real x_inv;
00316
00317 if (x.real_value == 0.0) return;
00318 x_inv = 1.0 / x.real_value;
00319
00320 if (fixed_axis) {
00321 ref1.set_weighted_gradient (-1.0 * x_inv * dist_v_ortho);
00322 main.set_weighted_gradient ( x_inv * dist_v_ortho);
00323 } else {
00324 if (b_no_PBC) {
00325 v13 = main.center_of_mass() - ref1.center_of_mass();
00326 } else {
00327 v13 = cvm::position_distance (ref1.center_of_mass(), main.center_of_mass());
00328 }
00329 A = (dist_v * axis) / axis_norm;
00330
00331 ref1.set_weighted_gradient ( (A - 1.0) * x_inv * dist_v_ortho);
00332 ref2.set_weighted_gradient ( -A * x_inv * dist_v_ortho);
00333 main.set_weighted_gradient ( 1.0 * x_inv * dist_v_ortho);
00334 }
00335 }
00336
00337 void colvar::distance_xy::calc_force_invgrads()
00338 {
00339 main.read_system_forces();
00340
00341 if (fixed_axis && !b_1site_force) {
00342 ref1.read_system_forces();
00343 ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho);
00344 } else {
00345 ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho;
00346 }
00347 }
00348
00349 void colvar::distance_xy::calc_Jacobian_derivative()
00350 {
00351 jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
00352 }
00353
00354 void colvar::distance_xy::apply_force (colvarvalue const &force)
00355 {
00356 if (!ref1.noforce)
00357 ref1.apply_colvar_force (force.real_value);
00358
00359 if (ref2.size() && !ref2.noforce)
00360 ref2.apply_colvar_force (force.real_value);
00361
00362 if (!main.noforce)
00363 main.apply_colvar_force (force.real_value);
00364 }
00365
00366
00367
00368 colvar::distance_dir::distance_dir (std::string const &conf)
00369 : distance (conf)
00370 {
00371 function_type = "distance_dir";
00372 x.type (colvarvalue::type_unitvector);
00373 }
00374
00375
00376 colvar::distance_dir::distance_dir()
00377 : distance()
00378 {
00379 function_type = "distance_dir";
00380 x.type (colvarvalue::type_unitvector);
00381 }
00382
00383
00384 void colvar::distance_dir::calc_value()
00385 {
00386 if (b_no_PBC) {
00387 dist_v = group2.center_of_mass() - group1.center_of_mass();
00388 } else {
00389 dist_v = cvm::position_distance (group1.center_of_mass(),
00390 group2.center_of_mass());
00391 }
00392 x.rvector_value = dist_v.unit();
00393 }
00394
00395
00396 void colvar::distance_dir::calc_gradients()
00397 {
00398
00399
00400
00401 }
00402
00403
00404 void colvar::distance_dir::apply_force (colvarvalue const &force)
00405 {
00406
00407 cvm::real const iprod = force.rvector_value * x.rvector_value;
00408 cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;
00409
00410 if (!group1.noforce)
00411 group1.apply_force (-1.0 * force_tang);
00412
00413 if (!group2.noforce)
00414 group2.apply_force ( force_tang);
00415 }
00416
00417
00418
00419 colvar::distance_inv::distance_inv (std::string const &conf)
00420 : distance (conf)
00421 {
00422 function_type = "distance_inv";
00423 get_keyval (conf, "exponent", exponent, 6);
00424 if (exponent%2) {
00425 cvm::fatal_error ("Error: odd exponent provided, can only use even ones.\n");
00426 }
00427 if (exponent <= 0) {
00428 cvm::fatal_error ("Error: negative or zero exponent provided.\n");
00429 }
00430
00431 for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
00432 for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
00433 if (ai1->id == ai2->id)
00434 cvm::fatal_error ("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n");
00435 }
00436 }
00437
00438 b_inverse_gradients = false;
00439 b_Jacobian_derivative = false;
00440 x.type (colvarvalue::type_scalar);
00441 }
00442
00443 colvar::distance_inv::distance_inv()
00444 {
00445 function_type = "distance_inv";
00446 exponent = 6;
00447 b_inverse_gradients = false;
00448 b_Jacobian_derivative = false;
00449 b_1site_force = false;
00450 x.type (colvarvalue::type_scalar);
00451 }
00452
00453 void colvar::distance_inv::calc_value()
00454 {
00455 x.real_value = 0.0;
00456 if (b_no_PBC) {
00457 for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
00458 for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
00459 cvm::rvector const dv = ai2->pos - ai1->pos;
00460 cvm::real const d2 = dv.norm2();
00461 cvm::real dinv = 1.0;
00462 for (int ne = 0; ne < exponent/2; ne++)
00463 dinv *= 1.0/d2;
00464 x.real_value += dinv;
00465 cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
00466 ai1->grad += -1.0 * dsumddv;
00467 ai2->grad += dsumddv;
00468 }
00469 }
00470 } else {
00471 for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
00472 for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
00473 cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
00474 cvm::real const d2 = dv.norm2();
00475 cvm::real dinv = 1.0;
00476 for (int ne = 0; ne < exponent/2; ne++)
00477 dinv *= 1.0/d2;
00478 x.real_value += dinv;
00479 cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
00480 ai1->grad += -1.0 * dsumddv;
00481 ai2->grad += dsumddv;
00482 }
00483 }
00484 }
00485
00486 x.real_value *= 1.0 / cvm::real (group1.size() * group2.size());
00487 x.real_value = std::pow (x.real_value, -1.0/(cvm::real (exponent)));
00488 }
00489
00490 void colvar::distance_inv::calc_gradients()
00491 {
00492 cvm::real const dxdsum = (-1.0/(cvm::real (exponent))) * std::pow (x.real_value, exponent+1) / cvm::real (group1.size() * group2.size());
00493 for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
00494 ai1->grad *= dxdsum;
00495 }
00496 for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
00497 ai2->grad *= dxdsum;
00498 }
00499 }
00500
00501 void colvar::distance_inv::apply_force (colvarvalue const &force)
00502 {
00503 if (!group1.noforce)
00504 group1.apply_colvar_force (force.real_value);
00505
00506 if (!group2.noforce)
00507 group2.apply_colvar_force (force.real_value);
00508 }
00509
00510
00511
00512 colvar::gyration::gyration (std::string const &conf)
00513 : cvc (conf)
00514 {
00515 function_type = "gyration";
00516 b_inverse_gradients = true;
00517 b_Jacobian_derivative = true;
00518 parse_group (conf, "atoms", atoms);
00519 atom_groups.push_back (&atoms);
00520
00521 if (atoms.b_user_defined_fit) {
00522 cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
00523 } else {
00524 atoms.b_center = true;
00525 atoms.ref_pos.assign (1, cvm::rvector (0.0, 0.0, 0.0));
00526 }
00527
00528 x.type (colvarvalue::type_scalar);
00529 }
00530
00531
00532 colvar::gyration::gyration()
00533 {
00534 function_type = "gyration";
00535 b_inverse_gradients = true;
00536 b_Jacobian_derivative = true;
00537 x.type (colvarvalue::type_scalar);
00538 }
00539
00540
00541 void colvar::gyration::calc_value()
00542 {
00543 x.real_value = 0.0;
00544 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00545 x.real_value += (ai->pos).norm2();
00546 }
00547 x.real_value = std::sqrt (x.real_value / cvm::real (atoms.size()));
00548 }
00549
00550
00551 void colvar::gyration::calc_gradients()
00552 {
00553 cvm::real const drdx = 1.0/(cvm::real (atoms.size()) * x.real_value);
00554 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00555 ai->grad = drdx * ai->pos;
00556 }
00557
00558 if (b_debug_gradients) {
00559 cvm::log ("Debugging gradients:\n");
00560 debug_gradients (atoms);
00561 }
00562 }
00563
00564
00565 void colvar::gyration::calc_force_invgrads()
00566 {
00567 atoms.read_system_forces();
00568
00569 cvm::real const dxdr = 1.0/x.real_value;
00570 ft.real_value = 0.0;
00571
00572 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00573 ft.real_value += dxdr * ai->pos * ai->system_force;
00574 }
00575 }
00576
00577
00578 void colvar::gyration::calc_Jacobian_derivative()
00579 {
00580 jd = x.real_value ? (3.0 * cvm::real (atoms.size()) - 4.0) / x.real_value : 0.0;
00581 }
00582
00583
00584 void colvar::gyration::apply_force (colvarvalue const &force)
00585 {
00586 if (!atoms.noforce)
00587 atoms.apply_colvar_force (force.real_value);
00588 }
00589
00590
00591
00592 colvar::inertia::inertia (std::string const &conf)
00593 : gyration (conf)
00594 {
00595 function_type = "inertia";
00596 b_inverse_gradients = false;
00597 b_Jacobian_derivative = false;
00598 x.type (colvarvalue::type_scalar);
00599 }
00600
00601
00602 colvar::inertia::inertia()
00603 {
00604 function_type = "inertia";
00605 x.type (colvarvalue::type_scalar);
00606 }
00607
00608
00609 void colvar::inertia::calc_value()
00610 {
00611 x.real_value = 0.0;
00612 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00613 x.real_value += (ai->pos).norm2();
00614 }
00615 }
00616
00617
00618 void colvar::inertia::calc_gradients()
00619 {
00620 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00621 ai->grad = 2.0 * ai->pos;
00622 }
00623
00624 if (b_debug_gradients) {
00625 cvm::log ("Debugging gradients:\n");
00626 debug_gradients (atoms);
00627 }
00628 }
00629
00630
00631 void colvar::inertia::apply_force (colvarvalue const &force)
00632 {
00633 if (!atoms.noforce)
00634 atoms.apply_colvar_force (force.real_value);
00635 }
00636
00637
00638 colvar::inertia_z::inertia_z (std::string const &conf)
00639 : inertia (conf)
00640 {
00641 function_type = "inertia_z";
00642 if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
00643 if (axis.norm2() == 0.0)
00644 cvm::fatal_error ("Axis vector is zero!");
00645 if (axis.norm2() != 1.0) {
00646 axis = axis.unit();
00647 cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
00648 }
00649 }
00650 x.type (colvarvalue::type_scalar);
00651 }
00652
00653
00654 colvar::inertia_z::inertia_z()
00655 {
00656 function_type = "inertia_z";
00657 x.type (colvarvalue::type_scalar);
00658 }
00659
00660
00661 void colvar::inertia_z::calc_value()
00662 {
00663 x.real_value = 0.0;
00664 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00665 cvm::real const iprod = ai->pos * axis;
00666 x.real_value += iprod * iprod;
00667 }
00668 }
00669
00670
00671 void colvar::inertia_z::calc_gradients()
00672 {
00673 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00674 ai->grad = 2.0 * (ai->pos * axis) * axis;
00675 }
00676
00677 if (b_debug_gradients) {
00678 cvm::log ("Debugging gradients:\n");
00679 debug_gradients (atoms);
00680 }
00681 }
00682
00683
00684 void colvar::inertia_z::apply_force (colvarvalue const &force)
00685 {
00686 if (!atoms.noforce)
00687 atoms.apply_colvar_force (force.real_value);
00688 }
00689
00690
00691
00692 colvar::rmsd::rmsd (std::string const &conf)
00693 : cvc (conf)
00694 {
00695 b_inverse_gradients = true;
00696 b_Jacobian_derivative = true;
00697 function_type = "rmsd";
00698 x.type (colvarvalue::type_scalar);
00699
00700 parse_group (conf, "atoms", atoms);
00701 atom_groups.push_back (&atoms);
00702
00703 if (atoms.b_dummy)
00704 cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom.");
00705
00706 if (atoms.ref_pos_group != NULL) {
00707 cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
00708 "Jacobian derivatives of the RMSD will not be calculated.\n");
00709 b_Jacobian_derivative = false;
00710 }
00711
00712
00713
00714
00715 if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
00716 cvm::log ("Using reference positions from configuration file to calculate the variable.\n");
00717 if (ref_pos.size() != atoms.size()) {
00718 cvm::fatal_error ("Error: the number of reference positions provided ("+
00719 cvm::to_str (ref_pos.size())+
00720 ") does not match the number of atoms of group \"atoms\" ("+
00721 cvm::to_str (atoms.size())+").\n");
00722 }
00723 }
00724 {
00725 std::string ref_pos_file;
00726 if (get_keyval (conf, "refPositionsFile", ref_pos_file, std::string (""))) {
00727
00728 if (ref_pos.size()) {
00729 cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and "
00730 "\"refPositions\" at the same time.\n");
00731 }
00732
00733 std::string ref_pos_col;
00734 double ref_pos_col_value;
00735
00736 if (get_keyval (conf, "refPositionsCol", ref_pos_col, std::string (""))) {
00737
00738 bool found = get_keyval (conf, "refPositionsColValue", ref_pos_col_value, 0.0);
00739 if (found && !ref_pos_col_value)
00740 cvm::fatal_error ("Error: refPositionsColValue, "
00741 "if provided, must be non-zero.\n");
00742 } else {
00743
00744 atoms.create_sorted_ids();
00745 }
00746
00747 ref_pos.resize (atoms.size());
00748 cvm::load_coords (ref_pos_file.c_str(), ref_pos, atoms.sorted_ids,
00749 ref_pos_col, ref_pos_col_value);
00750 }
00751 }
00752
00753 if (atoms.b_user_defined_fit) {
00754 cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
00755 } else {
00756
00757 cvm::log ("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: "
00758 "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
00759 atoms.b_center = true;
00760 atoms.b_rotate = true;
00761
00762
00763 atoms.ref_pos = ref_pos;
00764 atoms.center_ref_pos();
00765
00766
00767 atoms.rot.request_group1_gradients (atoms.size());
00768
00769
00770 atoms.rot.request_group2_gradients (atoms.size());
00771 }
00772 }
00773
00774
00775 void colvar::rmsd::calc_value()
00776 {
00777
00778
00779 x.real_value = 0.0;
00780 for (size_t ia = 0; ia < atoms.size(); ia++) {
00781 x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2();
00782 }
00783 x.real_value /= cvm::real (atoms.size());
00784 x.real_value = std::sqrt (x.real_value);
00785 }
00786
00787
00788 void colvar::rmsd::calc_gradients()
00789 {
00790 cvm::real const drmsddx2 = (x.real_value > 0.0) ?
00791 0.5 / (x.real_value * cvm::real (atoms.size())) :
00792 0.0;
00793
00794 for (size_t ia = 0; ia < atoms.size(); ia++) {
00795 atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - ref_pos[ia]));
00796 }
00797
00798 if (b_debug_gradients) {
00799 cvm::log ("Debugging gradients:\n");
00800 debug_gradients (atoms);
00801 }
00802 }
00803
00804
00805 void colvar::rmsd::apply_force (colvarvalue const &force)
00806 {
00807 if (!atoms.noforce)
00808 atoms.apply_colvar_force (force.real_value);
00809 }
00810
00811
00812 void colvar::rmsd::calc_force_invgrads()
00813 {
00814 atoms.read_system_forces();
00815 ft.real_value = 0.0;
00816
00817
00818
00819 for (size_t ia = 0; ia < atoms.size(); ia++) {
00820 ft.real_value += atoms[ia].grad * atoms[ia].system_force;
00821 }
00822 ft.real_value *= atoms.size();
00823 }
00824
00825
00826 void colvar::rmsd::calc_Jacobian_derivative()
00827 {
00828
00829 cvm::real divergence = 0.0;
00830
00831 if (atoms.b_rotate) {
00832
00833
00834 cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
00835
00836 cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
00837 for (size_t ia = 0; ia < atoms.size(); ia++) {
00838
00839
00840 cvm::vector1d< cvm::rvector, 4 > &dq = atoms.rot.dQ0_1[ia];
00841
00842 g11 = 2.0 * (atoms.rot.q)[1]*dq[1];
00843 g22 = 2.0 * (atoms.rot.q)[2]*dq[2];
00844 g33 = 2.0 * (atoms.rot.q)[3]*dq[3];
00845 g01 = (atoms.rot.q)[0]*dq[1] + (atoms.rot.q)[1]*dq[0];
00846 g02 = (atoms.rot.q)[0]*dq[2] + (atoms.rot.q)[2]*dq[0];
00847 g03 = (atoms.rot.q)[0]*dq[3] + (atoms.rot.q)[3]*dq[0];
00848 g12 = (atoms.rot.q)[1]*dq[2] + (atoms.rot.q)[2]*dq[1];
00849 g13 = (atoms.rot.q)[1]*dq[3] + (atoms.rot.q)[3]*dq[1];
00850 g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2];
00851
00852
00853 grad_rot_mat[0][0] = -2.0 * (g22 + g33);
00854 grad_rot_mat[1][0] = 2.0 * (g12 + g03);
00855 grad_rot_mat[2][0] = 2.0 * (g13 - g02);
00856 grad_rot_mat[0][1] = 2.0 * (g12 - g03);
00857 grad_rot_mat[1][1] = -2.0 * (g11 + g33);
00858 grad_rot_mat[2][1] = 2.0 * (g01 + g23);
00859 grad_rot_mat[0][2] = 2.0 * (g02 + g13);
00860 grad_rot_mat[1][2] = 2.0 * (g23 - g01);
00861 grad_rot_mat[2][2] = -2.0 * (g11 + g22);
00862
00863 cvm::atom_pos &y = ref_pos[ia];
00864
00865 for (size_t alpha = 0; alpha < 3; alpha++) {
00866 for (size_t beta = 0; beta < 3; beta++) {
00867 divergence += grad_rot_mat[beta][alpha][alpha] * y[beta];
00868
00869
00870
00871 }
00872 }
00873 }
00874 }
00875
00876 jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0;
00877 }
00878
00879
00880
00881
00882 colvar::eigenvector::eigenvector (std::string const &conf)
00883 : cvc (conf)
00884 {
00885 b_inverse_gradients = true;
00886 b_Jacobian_derivative = true;
00887 function_type = "eigenvector";
00888 x.type (colvarvalue::type_scalar);
00889
00890 parse_group (conf, "atoms", atoms);
00891 atom_groups.push_back (&atoms);
00892
00893
00894 {
00895 bool const b_inline = get_keyval (conf, "refPositions", ref_pos, ref_pos);
00896
00897 if (b_inline) {
00898 cvm::log ("Using reference positions from input file.\n");
00899 if (ref_pos.size() != atoms.size()) {
00900 cvm::fatal_error ("Error: reference positions do not "
00901 "match the number of requested atoms.\n");
00902 }
00903 }
00904
00905 std::string file_name;
00906 if (get_keyval (conf, "refPositionsFile", file_name)) {
00907
00908 if (b_inline)
00909 cvm::fatal_error ("Error: refPositions and refPositionsFile cannot be specified at the same time.\n");
00910
00911 std::string file_col;
00912 double file_col_value;
00913 if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) {
00914
00915 bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0);
00916 if (found && !file_col_value)
00917 cvm::fatal_error ("Error: refPositionsColValue, "
00918 "if provided, must be non-zero.\n");
00919 } else {
00920
00921 atoms.create_sorted_ids();
00922 }
00923
00924 ref_pos.resize (atoms.size());
00925 cvm::load_coords (file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value);
00926 }
00927 }
00928
00929
00930 cvm::rvector ref_pos_center (0.0, 0.0, 0.0);
00931 for (size_t i = 0; i < atoms.size(); i++) {
00932 ref_pos_center += ref_pos[i];
00933 }
00934 ref_pos_center *= 1.0 / atoms.size();
00935
00936 if (atoms.b_user_defined_fit) {
00937 cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
00938 } else {
00939
00940 cvm::log ("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: "
00941 "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
00942 atoms.b_center = true;
00943 atoms.b_rotate = true;
00944 atoms.ref_pos = ref_pos;
00945 atoms.center_ref_pos();
00946
00947
00948 atoms.rot.request_group1_gradients (atoms.size());
00949
00950
00951 atoms.rot.request_group2_gradients (atoms.size());
00952 }
00953
00954 {
00955 bool const b_inline = get_keyval (conf, "vector", eigenvec, eigenvec);
00956
00957 if (b_inline) {
00958 cvm::log ("Using vector components from input file.\n");
00959 if (eigenvec.size() != atoms.size()) {
00960 cvm::fatal_error ("Error: vector components do not "
00961 "match the number of requested atoms.\n");
00962 }
00963 }
00964
00965 std::string file_name;
00966 if (get_keyval (conf, "vectorFile", file_name)) {
00967
00968 if (b_inline)
00969 cvm::fatal_error ("Error: vector and vectorFile cannot be specified at the same time.\n");
00970
00971 std::string file_col;
00972 double file_col_value;
00973 if (get_keyval (conf, "vectorCol", file_col, std::string (""))) {
00974
00975 bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0);
00976 if (found && !file_col_value)
00977 cvm::fatal_error ("Error: vectorColValue, "
00978 "if provided, must be non-zero.\n");
00979 } else {
00980
00981 atoms.create_sorted_ids();
00982 }
00983
00984 eigenvec.resize (atoms.size());
00985 cvm::load_coords (file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value);
00986 }
00987 }
00988
00989 if (!ref_pos.size() || !eigenvec.size()) {
00990 cvm::fatal_error ("Error: both reference coordinates"
00991 "and eigenvector must be defined.\n");
00992 }
00993
00994 cvm::rvector eig_center (0.0, 0.0, 0.0);
00995 for (size_t i = 0; i < atoms.size(); i++) {
00996 eig_center += eigenvec[i];
00997 }
00998 eig_center *= 1.0 / atoms.size();
00999 cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n");
01000
01001 bool b_difference_vector = false;
01002 get_keyval (conf, "differenceVector", b_difference_vector, false);
01003
01004 if (b_difference_vector) {
01005
01006 if (atoms.b_center) {
01007
01008 for (size_t i = 0; i < atoms.size(); i++) {
01009 eigenvec[i] -= eig_center;
01010 ref_pos[i] -= ref_pos_center;
01011 }
01012 }
01013 if (atoms.b_rotate) {
01014 atoms.rot.calc_optimal_rotation (eigenvec, ref_pos);
01015 for (size_t i = 0; i < atoms.size(); i++) {
01016 eigenvec[i] = atoms.rot.rotate (eigenvec[i]);
01017 }
01018 }
01019 cvm::log ("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n");
01020 for (size_t i = 0; i < atoms.size(); i++) {
01021 eigenvec[i] -= ref_pos[i];
01022 }
01023 if (atoms.b_center) {
01024
01025 for (size_t i = 0; i < atoms.size(); i++) {
01026 ref_pos[i] += ref_pos_center;
01027 }
01028 }
01029
01030 } else {
01031 cvm::log ("Centering the provided vector to zero.\n");
01032 for (size_t i = 0; i < atoms.size(); i++) {
01033 eigenvec[i] -= eig_center;
01034 }
01035 }
01036
01037
01038
01039
01040 eigenvec_invnorm2 = 0.0;
01041 for (size_t i = 0; i < atoms.size(); i++) {
01042 eigenvec_invnorm2 += eigenvec[i].norm2();
01043 }
01044 eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
01045
01046 if (b_difference_vector) {
01047 cvm::log ("\"differenceVector\" is on: normalizing the vector.\n");
01048 for (size_t i = 0; i < atoms.size(); i++) {
01049 eigenvec[i] *= eigenvec_invnorm2;
01050 }
01051 } else {
01052 cvm::log ("The norm of the vector is |v| = "+cvm::to_str (eigenvec_invnorm2)+".\n");
01053 }
01054 }
01055
01056
01057 void colvar::eigenvector::calc_value()
01058 {
01059 x.real_value = 0.0;
01060 for (size_t i = 0; i < atoms.size(); i++) {
01061 x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i];
01062 }
01063 }
01064
01065
01066 void colvar::eigenvector::calc_gradients()
01067 {
01068 for (size_t ia = 0; ia < atoms.size(); ia++) {
01069 atoms[ia].grad = eigenvec[ia];
01070 }
01071
01072 if (b_debug_gradients) {
01073 cvm::log ("Debugging gradients:\n");
01074 debug_gradients (atoms);
01075 }
01076 }
01077
01078
01079 void colvar::eigenvector::apply_force (colvarvalue const &force)
01080 {
01081 if (!atoms.noforce)
01082 atoms.apply_colvar_force (force.real_value);
01083 }
01084
01085
01086 void colvar::eigenvector::calc_force_invgrads()
01087 {
01088 atoms.read_system_forces();
01089 ft.real_value = 0.0;
01090
01091 for (size_t ia = 0; ia < atoms.size(); ia++) {
01092 ft.real_value += eigenvec_invnorm2 * atoms[ia].grad *
01093 atoms[ia].system_force;
01094 }
01095 }
01096
01097
01098 void colvar::eigenvector::calc_Jacobian_derivative()
01099 {
01100
01101 cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
01102 cvm::quaternion &quat0 = atoms.rot.q;
01103
01104
01105 cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
01106
01107 cvm::atom_pos x_relative;
01108 cvm::real sum = 0.0;
01109
01110 for (size_t ia = 0; ia < atoms.size(); ia++) {
01111
01112
01113
01114
01115 cvm::vector1d< cvm::rvector, 4 > &dq_1 = atoms.rot.dQ0_1[ia];
01116
01117 g11 = 2.0 * quat0[1]*dq_1[1];
01118 g22 = 2.0 * quat0[2]*dq_1[2];
01119 g33 = 2.0 * quat0[3]*dq_1[3];
01120 g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0];
01121 g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0];
01122 g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0];
01123 g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1];
01124 g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1];
01125 g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2];
01126
01127
01128
01129 grad_rot_mat[0][0] = -2.0 * (g22 + g33);
01130 grad_rot_mat[0][1] = 2.0 * (g12 + g03);
01131 grad_rot_mat[0][2] = 2.0 * (g13 - g02);
01132 grad_rot_mat[1][0] = 2.0 * (g12 - g03);
01133 grad_rot_mat[1][1] = -2.0 * (g11 + g33);
01134 grad_rot_mat[1][2] = 2.0 * (g01 + g23);
01135 grad_rot_mat[2][0] = 2.0 * (g02 + g13);
01136 grad_rot_mat[2][1] = 2.0 * (g23 - g01);
01137 grad_rot_mat[2][2] = -2.0 * (g11 + g22);
01138
01139 for (size_t i = 0; i < 3; i++) {
01140 for (size_t j = 0; j < 3; j++) {
01141 sum += grad_rot_mat[i][j][i] * eigenvec[ia][j];
01142 }
01143 }
01144 }
01145
01146 jd.real_value = sum * std::sqrt (eigenvec_invnorm2);
01147 }
01148
01149
01150
01151