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
00013
00014
00015 colvar::distance::distance (std::string const &conf)
00016 : cvc (conf)
00017 {
00018 function_type = "distance";
00019 b_inverse_gradients = true;
00020 b_Jacobian_derivative = true;
00021 if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
00022 cvm::log ("Computing system force on group 1 only");
00023 }
00024 parse_group (conf, "group1", group1);
00025 parse_group (conf, "group2", group2);
00026 x.type (colvarvalue::type_scalar);
00027 }
00028
00029 colvar::distance::distance()
00030 : cvc ()
00031 {
00032 function_type = "distance";
00033 b_inverse_gradients = true;
00034 b_Jacobian_derivative = true;
00035 b_1site_force = false;
00036 x.type (colvarvalue::type_scalar);
00037 }
00038
00039 void colvar::distance::calc_value()
00040 {
00041 group1.reset_atoms_data();
00042 group2.reset_atoms_data();
00043
00044 group1.read_positions();
00045 group2.read_positions();
00046
00047 dist_v = cvm::position_distance (group1.center_of_mass(),
00048 group2.center_of_mass());
00049 x.real_value = dist_v.norm();
00050 }
00051
00052 void colvar::distance::calc_gradients()
00053 {
00054 cvm::rvector const u = dist_v.unit();
00055 group1.set_weighted_gradient (-1.0 * u);
00056 group2.set_weighted_gradient ( u);
00057 }
00058
00059 void colvar::distance::calc_force_invgrads()
00060 {
00061 group1.read_system_forces();
00062 if ( b_1site_force ) {
00063 ft.real_value = -1.0 * (group1.system_force() * dist_v.unit());
00064 } else {
00065 group2.read_system_forces();
00066 ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit());
00067 }
00068 }
00069
00070 void colvar::distance::calc_Jacobian_derivative()
00071 {
00072 jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;
00073 }
00074
00075 void colvar::distance::apply_force (colvarvalue const &force)
00076 {
00077 if (!group1.noforce)
00078 group1.apply_colvar_force (force);
00079
00080 if (!group2.noforce)
00081 group2.apply_colvar_force (force);
00082 }
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 colvar::distance_z::distance_z (std::string const &conf)
00132 : cvc (conf)
00133 {
00134 function_type = "distance_z";
00135 b_inverse_gradients = true;
00136 b_Jacobian_derivative = true;
00137 x.type (colvarvalue::type_scalar);
00138
00139 parse_group (conf, "main", main);
00140 parse_group (conf, "ref", ref1);
00141
00142 parse_group (conf, "ref2", ref2, true);
00143
00144 if (ref2.size()) {
00145 cvm::log ("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
00146 fixed_axis = false;
00147 if (key_lookup (conf, "axis"))
00148 cvm::log ("Warning: explicit axis definition will be ignored!");
00149 } else {
00150 if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
00151 if (axis.norm2() == 0.0)
00152 cvm::fatal_error ("Axis vector is zero!");
00153 axis = axis.unit();
00154 }
00155 fixed_axis = true;
00156 }
00157
00158 if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
00159 cvm::log ("Computing system force on group \"main\" only");
00160 }
00161 }
00162
00163 colvar::distance_z::distance_z()
00164 {
00165 function_type = "distance_z";
00166 b_inverse_gradients = true;
00167 b_Jacobian_derivative = true;
00168 x.type (colvarvalue::type_scalar);
00169 }
00170
00171 void colvar::distance_z::calc_value()
00172 {
00173 main.reset_atoms_data();
00174 ref1.reset_atoms_data();
00175
00176 main.read_positions();
00177 ref1.read_positions();
00178
00179 if (fixed_axis) {
00180 dist_v = cvm::position_distance (ref1.center_of_mass(),
00181 main.center_of_mass());
00182 } else {
00183 ref2.reset_atoms_data();
00184 ref2.read_positions();
00185
00186 dist_v = cvm::position_distance (0.5 * (ref1.center_of_mass() + ref2.center_of_mass()),
00187 main.center_of_mass());
00188 axis = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
00189 axis_norm = axis.norm();
00190 axis = axis.unit();
00191 }
00192 x.real_value = axis * dist_v;
00193 }
00194
00195 void colvar::distance_z::calc_gradients()
00196 {
00197 if (fixed_axis) {
00198 ref1.set_weighted_gradient (-1.0 * axis);
00199 main.set_weighted_gradient ( axis);
00200 } else {
00201 ref1.set_weighted_gradient ( 1.0 / axis_norm * (
00202 cvm::position_distance (ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis ));
00203 ref2.set_weighted_gradient ( 1.0 / axis_norm * (
00204 cvm::position_distance (main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis ));
00205 main.set_weighted_gradient ( axis );
00206 }
00207 }
00208
00209 void colvar::distance_z::calc_force_invgrads()
00210 {
00211 main.read_system_forces();
00212
00213 if (fixed_axis && !b_1site_force) {
00214 ref1.read_system_forces();
00215 ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis);
00216 } else {
00217 ft.real_value = main.system_force() * axis;
00218 }
00219 }
00220
00221 void colvar::distance_z::calc_Jacobian_derivative()
00222 {
00223 jd.real_value = 0.0;
00224 }
00225
00226 void colvar::distance_z::apply_force (colvarvalue const &force)
00227 {
00228 if (!ref1.noforce)
00229 ref1.apply_colvar_force (force.real_value);
00230
00231 if (ref2.size() && !ref2.noforce)
00232 ref2.apply_colvar_force (force.real_value);
00233
00234 if (!main.noforce)
00235 main.apply_colvar_force (force.real_value);
00236 }
00237
00238
00239
00240 colvar::distance_xy::distance_xy (std::string const &conf)
00241 : distance_z (conf)
00242 {
00243 function_type = "distance_xy";
00244 b_inverse_gradients = true;
00245 b_Jacobian_derivative = true;
00246 x.type (colvarvalue::type_scalar);
00247 }
00248
00249 colvar::distance_xy::distance_xy()
00250 : distance_z()
00251 {
00252 function_type = "distance_xy";
00253 b_inverse_gradients = true;
00254 b_Jacobian_derivative = true;
00255 x.type (colvarvalue::type_scalar);
00256 }
00257
00258 void colvar::distance_xy::calc_value()
00259 {
00260 ref1.reset_atoms_data();
00261 main.reset_atoms_data();
00262
00263 ref1.read_positions();
00264 main.read_positions();
00265
00266 dist_v = cvm::position_distance (ref1.center_of_mass(),
00267 main.center_of_mass());
00268 if (!fixed_axis) {
00269 ref2.reset_atoms_data();
00270 ref2.read_positions();
00271
00272 v12 = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
00273 axis_norm = v12.norm();
00274 axis = v12.unit();
00275 }
00276
00277 dist_v_ortho = dist_v - (dist_v * axis) * axis;
00278 x.real_value = dist_v_ortho.norm();
00279 }
00280
00281 void colvar::distance_xy::calc_gradients()
00282 {
00283
00284
00285 cvm::real A;
00286 cvm::real x_inv;
00287
00288 if (x.real_value == 0.0) return;
00289 x_inv = 1.0 / x.real_value;
00290
00291 if (fixed_axis) {
00292 ref1.set_weighted_gradient (-1.0 * x_inv * dist_v_ortho);
00293 main.set_weighted_gradient ( x_inv * dist_v_ortho);
00294 } else {
00295 v13 = cvm::position_distance (ref1.center_of_mass(), main.center_of_mass());
00296 A = (dist_v * axis) / axis_norm;
00297
00298 ref1.set_weighted_gradient ( (A - 1.0) * x_inv * dist_v_ortho);
00299 ref2.set_weighted_gradient ( -A * x_inv * dist_v_ortho);
00300 main.set_weighted_gradient ( 1.0 * x_inv * dist_v_ortho);
00301 }
00302 }
00303
00304 void colvar::distance_xy::calc_force_invgrads()
00305 {
00306 main.read_system_forces();
00307
00308 if (fixed_axis && !b_1site_force) {
00309 ref1.read_system_forces();
00310 ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho);
00311 } else {
00312 ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho;
00313 }
00314 }
00315
00316 void colvar::distance_xy::calc_Jacobian_derivative()
00317 {
00318 jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
00319 }
00320
00321 void colvar::distance_xy::apply_force (colvarvalue const &force)
00322 {
00323 if (!ref1.noforce)
00324 ref1.apply_colvar_force (force.real_value);
00325
00326 if (ref2.size() && !ref2.noforce)
00327 ref2.apply_colvar_force (force.real_value);
00328
00329 if (!main.noforce)
00330 main.apply_colvar_force (force.real_value);
00331 }
00332
00333
00334
00335 colvar::min_distance::min_distance (std::string const &conf)
00336 : distance (conf)
00337 {
00338 function_type = "min_distance";
00339 x.type (colvarvalue::type_scalar);
00340
00341 get_keyval (conf, "smoothing", smoothing, (1.0 * cvm::unit_angstrom()));
00342 }
00343
00344 colvar::min_distance::min_distance()
00345 : distance()
00346 {
00347 function_type = "min_distance";
00348 x.type (colvarvalue::type_scalar);
00349 }
00350
00351 void colvar::min_distance::calc_value()
00352 {
00353 group1.reset_atoms_data();
00354 group2.reset_atoms_data();
00355
00356 group1.read_positions();
00357 group2.read_positions();
00358
00359 x.real_value = 0.0;
00360
00361 bool zero_dist = false;
00362
00363 for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
00364 for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
00365 cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
00366 cvm::real const d = dv.norm();
00367 if (d > 0.0)
00368 x.real_value += ::exp (smoothing / d);
00369 else
00370 zero_dist = true;
00371 }
00372 }
00373
00374 x.real_value = zero_dist ? 0.0 : smoothing/(::log (x.real_value));
00375 }
00376
00377 void colvar::min_distance::calc_gradients()
00378 {
00379 if (x.real_value > 0.0) {
00380 cvm::real const sum = ::exp (smoothing/x.real_value);
00381 cvm::real const dxdsum = -1.0 *
00382 (x.real_value/smoothing) * (x.real_value/smoothing) *
00383 (1.0 / sum);
00384
00385 for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
00386 for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
00387 cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
00388 cvm::real const d = dv.norm();
00389 if (d > 0.0) {
00390 cvm::rvector const dvu = dv / dv.norm();
00391 ai1->grad += dxdsum * ::exp (smoothing / d) *
00392 smoothing * (-1.0/(d*d)) * (-1.0) * dvu;
00393 ai2->grad += dxdsum * ::exp (smoothing / d) *
00394 smoothing * (-1.0/(d*d)) * dvu;
00395 }
00396 }
00397 }
00398 }
00399 }
00400
00401 void colvar::min_distance::apply_force (colvarvalue const &force)
00402 {
00403 if (!group1.noforce)
00404 group1.apply_colvar_force (force.real_value);
00405
00406 if (!group2.noforce)
00407 group2.apply_colvar_force (force.real_value);
00408 }
00409
00410
00411
00412 colvar::distance_dir::distance_dir (std::string const &conf)
00413 : distance (conf)
00414 {
00415 function_type = "distance_dir";
00416 x.type (colvarvalue::type_unitvector);
00417 }
00418
00419
00420 colvar::distance_dir::distance_dir()
00421 : distance()
00422 {
00423 function_type = "distance_dir";
00424 x.type (colvarvalue::type_unitvector);
00425 }
00426
00427
00428 void colvar::distance_dir::calc_value()
00429 {
00430 group1.reset_atoms_data();
00431 group2.reset_atoms_data();
00432
00433 group1.read_positions();
00434 group2.read_positions();
00435
00436 dist_v = cvm::position_distance (group1.center_of_mass(),
00437 group2.center_of_mass());
00438 x.rvector_value = dist_v.unit();
00439 }
00440
00441
00442 void colvar::distance_dir::calc_gradients()
00443 {
00444
00445 }
00446
00447
00448 void colvar::distance_dir::apply_force (colvarvalue const &force)
00449 {
00450
00451 cvm::real const iprod = force.rvector_value * x.rvector_value;
00452 cvm::rvector const force_tang = (1.0 - iprod) * force.rvector_value;
00453
00454 if (!group1.noforce)
00455 group1.apply_force (-1.0 * force_tang);
00456
00457 if (!group2.noforce)
00458 group2.apply_force ( force_tang);
00459 }
00460
00461
00462
00463
00464 colvar::gyration::gyration (std::string const &conf)
00465 : cvc (conf)
00466 {
00467 function_type = "gyration";
00468 b_inverse_gradients = true;
00469 b_Jacobian_derivative = true;
00470 parse_group (conf, "atoms", atoms);
00471 x.type (colvarvalue::type_scalar);
00472 }
00473
00474
00475 colvar::gyration::gyration()
00476 {
00477 function_type = "gyration";
00478 b_inverse_gradients = true;
00479 b_Jacobian_derivative = true;
00480 x.type (colvarvalue::type_scalar);
00481 }
00482
00483
00484 void colvar::gyration::calc_value()
00485 {
00486 atoms.reset_atoms_data();
00487 atoms.read_positions();
00488 atoms.apply_translation ((-1.0) * atoms.center_of_geometry());
00489
00490 x.real_value = 0.0;
00491 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00492 x.real_value += (ai->mass/atoms.total_mass) * (ai->pos).norm2();
00493 }
00494 x.real_value = ::sqrt (x.real_value);
00495 }
00496
00497
00498 void colvar::gyration::calc_gradients()
00499 {
00500 cvm::real const drdx = 1.0/(cvm::real (atoms.size()) * x.real_value);
00501 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00502 ai->grad = drdx * ai->pos;
00503 }
00504 }
00505
00506
00507 void colvar::gyration::calc_force_invgrads()
00508 {
00509 atoms.read_system_forces();
00510
00511 cvm::real const dxdr = 1.0/x.real_value;
00512 ft.real_value = 0.0;
00513
00514 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00515 ft.real_value += dxdr * ai->pos * ai->system_force;
00516 }
00517 }
00518
00519
00520 void colvar::gyration::calc_Jacobian_derivative()
00521 {
00522 jd = x.real_value ? (3.0 * cvm::real (atoms.size()) - 4.0) / x.real_value : 0.0;
00523 }
00524
00525
00526 void colvar::gyration::apply_force (colvarvalue const &force)
00527 {
00528 if (!atoms.noforce)
00529 atoms.apply_colvar_force (force.real_value);
00530 }
00531
00532
00533
00534 colvar::rmsd::rmsd (std::string const &conf)
00535 : orientation (conf)
00536 {
00537
00538
00539 b_inverse_gradients = true;
00540 b_Jacobian_derivative = true;
00541 function_type = "rmsd";
00542 x.type (colvarvalue::type_scalar);
00543
00544 ref_pos_sum2 = 0.0;
00545 for (size_t i = 0; i < ref_pos.size(); i++) {
00546 ref_pos_sum2 += ref_pos[i].norm2();
00547 }
00548 }
00549
00550
00551 void colvar::rmsd::calc_value()
00552 {
00553 atoms.reset_atoms_data();
00554 atoms.read_positions();
00555
00556 atoms_cog = atoms.center_of_geometry();
00557 rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
00558
00559 cvm::real group_pos_sum2 = 0.0;
00560 for (size_t i = 0; i < atoms.size(); i++) {
00561 group_pos_sum2 += (atoms[i].pos - atoms_cog).norm2();
00562 }
00563
00564
00565 cvm::real const MSD = 1.0/(cvm::real (atoms.size())) *
00566 ( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda );
00567
00568 x.real_value = (MSD > 0.0) ? ::sqrt (MSD) : 0.0;
00569 }
00570
00571
00572 void colvar::rmsd::calc_gradients()
00573 {
00574 cvm::real const drmsddx2 = (x.real_value > 0.0) ?
00575 0.5 / (x.real_value * cvm::real (atoms.size())) :
00576 0.0;
00577
00578 for (size_t ia = 0; ia < atoms.size(); ia++) {
00579 atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - atoms_cog -
00580 rot.q.rotate (ref_pos[ia])));
00581 }
00582 }
00583
00584
00585 void colvar::rmsd::apply_force (colvarvalue const &force)
00586 {
00587 if (!atoms.noforce)
00588 atoms.apply_colvar_force (force.real_value);
00589 }
00590
00591
00592 void colvar::rmsd::calc_force_invgrads()
00593 {
00594 atoms.read_system_forces();
00595 ft.real_value = 0.0;
00596
00597
00598
00599 for (size_t ia = 0; ia < atoms.size(); ia++) {
00600 ft.real_value += atoms[ia].grad * atoms[ia].system_force;
00601 }
00602 ft.real_value *= atoms.size();
00603 }
00604
00605
00606 void colvar::rmsd::calc_Jacobian_derivative()
00607 {
00608
00609 cvm::real divergence = 0.0;
00610
00611
00612 cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
00613
00614
00615 cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
00616
00617 for (size_t ia = 0; ia < atoms.size(); ia++) {
00618
00619
00620 cvm::vector1d< cvm::rvector, 4 > &dq = rot.dQ0_2[ia];
00621
00622 g11 = 2.0 * (rot.q)[1]*dq[1];
00623 g22 = 2.0 * (rot.q)[2]*dq[2];
00624 g33 = 2.0 * (rot.q)[3]*dq[3];
00625 g01 = (rot.q)[0]*dq[1] + (rot.q)[1]*dq[0];
00626 g02 = (rot.q)[0]*dq[2] + (rot.q)[2]*dq[0];
00627 g03 = (rot.q)[0]*dq[3] + (rot.q)[3]*dq[0];
00628 g12 = (rot.q)[1]*dq[2] + (rot.q)[2]*dq[1];
00629 g13 = (rot.q)[1]*dq[3] + (rot.q)[3]*dq[1];
00630 g23 = (rot.q)[2]*dq[3] + (rot.q)[3]*dq[2];
00631
00632
00633 grad_rot_mat[0][0] = -2.0 * (g22 + g33);
00634 grad_rot_mat[1][0] = 2.0 * (g12 + g03);
00635 grad_rot_mat[2][0] = 2.0 * (g13 - g02);
00636 grad_rot_mat[0][1] = 2.0 * (g12 - g03);
00637 grad_rot_mat[1][1] = -2.0 * (g11 + g33);
00638 grad_rot_mat[2][1] = 2.0 * (g01 + g23);
00639 grad_rot_mat[0][2] = 2.0 * (g02 + g13);
00640 grad_rot_mat[1][2] = 2.0 * (g23 - g01);
00641 grad_rot_mat[2][2] = -2.0 * (g11 + g22);
00642
00643 cvm::atom_pos &y = ref_pos[ia];
00644
00645 for (size_t i = 0; i < 3; i++) {
00646 for (size_t j = 0; j < 3; j++) {
00647 divergence += grad_rot_mat[i][j][i] * y[j];
00648 }
00649 }
00650 }
00651
00652 jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0;
00653 }
00654
00655
00656
00657 colvar::logmsd::logmsd (std::string const &conf)
00658 : orientation (conf)
00659 {
00660 b_inverse_gradients = true;
00661 b_Jacobian_derivative = true;
00662 function_type = "logmsd";
00663 x.type (colvarvalue::type_scalar);
00664
00665 ref_pos_sum2 = 0.0;
00666 for (size_t i = 0; i < ref_pos.size(); i++) {
00667 ref_pos_sum2 += ref_pos[i].norm2();
00668 }
00669 }
00670
00671
00672 void colvar::logmsd::calc_value()
00673 {
00674 atoms.reset_atoms_data();
00675 atoms.read_positions();
00676
00677 if (cvm::debug())
00678 cvm::log ("colvar::logmsd: current com: "+
00679 cvm::to_str (atoms.center_of_mass())+"\n");
00680
00681 atoms_cog = atoms.center_of_geometry();
00682 rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
00683
00684 cvm::real group_pos_sum2 = 0.0;
00685 for (size_t i = 0; i < atoms.size(); i++) {
00686 group_pos_sum2 += (atoms[i].pos-atoms_cog).norm2();
00687 }
00688
00689
00690 MSD = 1.0/(cvm::real (atoms.size())) *
00691 ( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda );
00692
00693 x.real_value = (MSD > 0.0) ? ::log(MSD) : 0.0;
00694 }
00695
00696
00697 void colvar::logmsd::calc_gradients()
00698 {
00699 cvm::real fact = (MSD > 0.0) ? 2.0/(cvm::real (atoms.size()) * MSD) : 0.0;
00700
00701 for (size_t ia = 0; ia < atoms.size(); ia++) {
00702 atoms[ia].grad = fact * (atoms[ia].pos - atoms_cog - rot.dL0_2[ia]);
00703 }
00704 }
00705
00706
00707 void colvar::logmsd::apply_force (colvarvalue const &force)
00708 {
00709 if (!atoms.noforce)
00710 atoms.apply_colvar_force (force.real_value);
00711 }
00712
00713
00714 void colvar::logmsd::calc_force_invgrads()
00715 {
00716 atoms.read_system_forces();
00717 ft.real_value = 0.0;
00718
00719
00720
00721 for (size_t ia = 0; ia < atoms.size(); ia++) {
00722 ft.real_value += atoms[ia].grad * atoms[ia].system_force;
00723 }
00724 ft.real_value *= atoms.size() * MSD / 4.0;
00725 }
00726
00727
00728 void colvar::logmsd::calc_Jacobian_derivative()
00729 {
00730
00731 cvm::real divergence = 0.0;
00732
00733
00734 cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
00735
00736
00737 cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
00738
00739 for (size_t ia = 0; ia < atoms.size(); ia++) {
00740
00741
00742 cvm::vector1d< cvm::rvector, 4 > &dq = rot.dQ0_2[ia];
00743
00744 g11 = 2.0 * (rot.q)[1]*dq[1];
00745 g22 = 2.0 * (rot.q)[2]*dq[2];
00746 g33 = 2.0 * (rot.q)[3]*dq[3];
00747 g01 = (rot.q)[0]*dq[1] + (rot.q)[1]*dq[0];
00748 g02 = (rot.q)[0]*dq[2] + (rot.q)[2]*dq[0];
00749 g03 = (rot.q)[0]*dq[3] + (rot.q)[3]*dq[0];
00750 g12 = (rot.q)[1]*dq[2] + (rot.q)[2]*dq[1];
00751 g13 = (rot.q)[1]*dq[3] + (rot.q)[3]*dq[1];
00752 g23 = (rot.q)[2]*dq[3] + (rot.q)[3]*dq[2];
00753
00754
00755
00756 grad_rot_mat[0][0] = -2.0 * (g22 + g33);
00757 grad_rot_mat[1][0] = 2.0 * (g12 + g03);
00758 grad_rot_mat[2][0] = 2.0 * (g13 - g02);
00759 grad_rot_mat[0][1] = 2.0 * (g12 - g03);
00760 grad_rot_mat[1][1] = -2.0 * (g11 + g33);
00761 grad_rot_mat[2][1] = 2.0 * (g01 + g23);
00762 grad_rot_mat[0][2] = 2.0 * (g02 + g13);
00763 grad_rot_mat[1][2] = 2.0 * (g23 - g01);
00764 grad_rot_mat[2][2] = -2.0 * (g11 + g22);
00765
00766 cvm::atom_pos &y = ref_pos[ia];
00767
00768 for (size_t i = 0; i < 3; i++) {
00769 for (size_t j = 0; j < 3; j++) {
00770 divergence += grad_rot_mat[i][j][i] * y[j];
00771 }
00772 }
00773 }
00774
00775 jd.real_value = (3.0 * atoms.size() - 3.0 - divergence) / 2.0;
00776 }
00777
00778
00779
00780 colvar::eigenvector::eigenvector (std::string const &conf)
00781 : cvc (conf)
00782 {
00783 b_inverse_gradients = true;
00784 b_Jacobian_derivative = true;
00785 function_type = "eigenvector";
00786 x.type (colvarvalue::type_scalar);
00787
00788 parse_group (conf, "atoms", atoms);
00789
00790 if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
00791 cvm::log ("Using reference positions from input file.\n");
00792 if (ref_pos.size() != atoms.size()) {
00793 cvm::fatal_error ("Error: reference positions do not "
00794 "match the number of atom indexes.\n");
00795 }
00796 }
00797
00798 {
00799 std::string file_name;
00800 if (get_keyval (conf, "refPositionsFile", file_name)) {
00801
00802 std::string file_col;
00803 get_keyval (conf, "refPositionsCol", file_col, std::string ("O"));
00804
00805 double file_col_value;
00806 bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0);
00807 if (found && !file_col_value)
00808 cvm::fatal_error ("Error: refPositionsColValue, "
00809 "if provided, must be non-zero.\n");
00810
00811 ref_pos.resize (atoms.size());
00812 cvm::load_coords (file_name.c_str(), ref_pos, file_col, file_col_value);
00813 }
00814 }
00815
00816
00817 if (get_keyval (conf, "vector", eigenvec, eigenvec)) {
00818 cvm::log ("Using reference positions from input file.\n");
00819 if (eigenvec.size() != atoms.size()) {
00820 cvm::fatal_error ("Error: reference positions do not "
00821 "match the number of atom indexes.\n");
00822 }
00823 }
00824
00825 {
00826 std::string file_name;
00827 if (get_keyval (conf, "vectorFile", file_name)) {
00828
00829 std::string file_col;
00830 get_keyval (conf, "vectorCol", file_col, std::string ("O"));
00831
00832 double file_col_value;
00833 bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0);
00834 if (found && !file_col_value)
00835 cvm::fatal_error ("Error: eigenvectorColValue, "
00836 "if provided, must be non-zero.\n");
00837
00838 eigenvec.resize (atoms.size());
00839 cvm::load_coords (file_name.c_str(), eigenvec, file_col, file_col_value);
00840 }
00841 }
00842
00843 if (!ref_pos.size() || !eigenvec.size()) {
00844 cvm::fatal_error ("Error: must define both reference "
00845 "coordinates and eigenvector.\n");
00846 }
00847
00848 cvm::rvector center (0.0, 0.0, 0.0);
00849 eigenvec_invnorm2 = 0.0;
00850
00851 for (size_t i = 0; i < atoms.size(); i++) {
00852 center += eigenvec[i];
00853 }
00854
00855 cvm::log ("Subtracting sum of eigenvector components: " + cvm::to_str (center) + "\n");
00856
00857 for (size_t i = 0; i < atoms.size(); i++) {
00858 eigenvec[i] = eigenvec[i] - center;
00859 eigenvec_invnorm2 += eigenvec[i].norm2();
00860 }
00861 eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
00862
00863
00864
00865 atoms.rot.request_group1_gradients(atoms.size());
00866 atoms.rot.request_group2_gradients(atoms.size());
00867 }
00868
00869
00870 void colvar::eigenvector::calc_value()
00871 {
00872 atoms.reset_atoms_data();
00873 atoms.read_positions();
00874
00875 x.real_value = 0.0;
00876 for (size_t i = 0; i < atoms.size(); i++) {
00877 x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i];
00878 }
00879 }
00880
00881
00882 void colvar::eigenvector::calc_gradients()
00883 {
00884 for (size_t i = 0; i < atoms.size(); i++) {
00885 atoms[i].grad = eigenvec[i];
00886 }
00887
00888
00889
00890 }
00891
00892
00893 void colvar::eigenvector::apply_force (colvarvalue const &force)
00894 {
00895 if (!atoms.noforce)
00896 atoms.apply_colvar_force (force.real_value);
00897 }
00898
00899
00900 void colvar::eigenvector::calc_force_invgrads()
00901 {
00902 atoms.read_system_forces();
00903 ft.real_value = 0.0;
00904
00905 for (size_t ia = 0; ia < atoms.size(); ia++) {
00906 ft.real_value += eigenvec_invnorm2 * atoms[ia].grad *
00907 atoms[ia].system_force;
00908 }
00909 }
00910
00911
00912 void colvar::eigenvector::calc_Jacobian_derivative()
00913 {
00914
00915 cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
00916 cvm::quaternion &quat0 = atoms.rot.q;
00917
00918
00919 cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
00920
00921 cvm::atom_pos x_relative;
00922 cvm::real sum = 0.0;
00923
00924 for (size_t ia = 0; ia < atoms.size(); ia++) {
00925
00926
00927
00928
00929 cvm::vector1d< cvm::rvector, 4 > &dq_1 = atoms.rot.dQ0_1[ia];
00930
00931 g11 = 2.0 * quat0[1]*dq_1[1];
00932 g22 = 2.0 * quat0[2]*dq_1[2];
00933 g33 = 2.0 * quat0[3]*dq_1[3];
00934 g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0];
00935 g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0];
00936 g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0];
00937 g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1];
00938 g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1];
00939 g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2];
00940
00941
00942
00943 grad_rot_mat[0][0] = -2.0 * (g22 + g33);
00944 grad_rot_mat[0][1] = 2.0 * (g12 + g03);
00945 grad_rot_mat[0][2] = 2.0 * (g13 - g02);
00946 grad_rot_mat[1][0] = 2.0 * (g12 - g03);
00947 grad_rot_mat[1][1] = -2.0 * (g11 + g33);
00948 grad_rot_mat[1][2] = 2.0 * (g01 + g23);
00949 grad_rot_mat[2][0] = 2.0 * (g02 + g13);
00950 grad_rot_mat[2][1] = 2.0 * (g23 - g01);
00951 grad_rot_mat[2][2] = -2.0 * (g11 + g22);
00952
00953 for (size_t i = 0; i < 3; i++) {
00954 for (size_t j = 0; j < 3; j++) {
00955 sum += grad_rot_mat[i][j][i] * eigenvec[ia][j];
00956 }
00957 }
00958 }
00959
00960 jd.real_value = sum * sqrt (eigenvec_invnorm2);
00961 }