version 1.15 | version 1.16 |
---|
| |
/// -*- c++ -*- | // -*- c++ -*- |
| |
#include <cmath> | #include <cmath> |
| |
| |
| |
| |
| |
// "twogroup" flag defaults to true; set to false by selfCoordNum | colvar::distance::distance(std::string const &conf) |
// (only distance-derived component based on only one group) | |
| |
colvar::distance::distance(std::string const &conf, bool twogroups) | |
: cvc(conf) | : cvc(conf) |
{ | { |
function_type = "distance"; | function_type = "distance"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
| provide(f_cvc_com_based); |
| |
| group1 = parse_group(conf, "group1"); |
| group2 = parse_group(conf, "group2"); |
| |
if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { | if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { |
cvm::log("Computing distance using absolute positions (not minimal-image)"); | cvm::log("Computing distance using absolute positions (not minimal-image)"); |
} | } |
if (twogroups && get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { | |
cvm::log("Computing system force on group 1 only"); | init_total_force_params(conf); |
} | |
parse_group(conf, "group1", group1); | |
atom_groups.push_back(&group1); | |
if (twogroups) { | |
parse_group(conf, "group2", group2); | |
atom_groups.push_back(&group2); | |
} | |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
: cvc() | : cvc() |
{ | { |
function_type = "distance"; | function_type = "distance"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
b_1site_force = false; | provide(f_cvc_com_based); |
b_no_PBC = false; | b_no_PBC = false; |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
void colvar::distance::calc_value() | void colvar::distance::calc_value() |
{ | { |
if (b_no_PBC) { | if (b_no_PBC) { |
dist_v = group2.center_of_mass() - group1.center_of_mass(); | dist_v = group2->center_of_mass() - group1->center_of_mass(); |
} else { | } else { |
dist_v = cvm::position_distance(group1.center_of_mass(), | dist_v = cvm::position_distance(group1->center_of_mass(), |
group2.center_of_mass()); | group2->center_of_mass()); |
} | } |
x.real_value = dist_v.norm(); | x.real_value = dist_v.norm(); |
} | } |
| |
| |
void colvar::distance::calc_gradients() | void colvar::distance::calc_gradients() |
{ | { |
cvm::rvector const u = dist_v.unit(); | cvm::rvector const u = dist_v.unit(); |
group1.set_weighted_gradient(-1.0 * u); | group1->set_weighted_gradient(-1.0 * u); |
group2.set_weighted_gradient( u); | group2->set_weighted_gradient( u); |
} | } |
| |
| |
void colvar::distance::calc_force_invgrads() | void colvar::distance::calc_force_invgrads() |
{ | { |
group1.read_system_forces(); | group1->read_total_forces(); |
if ( b_1site_force ) { | if (is_enabled(f_cvc_one_site_total_force)) { |
ft.real_value = -1.0 * (group1.system_force() * dist_v.unit()); | ft.real_value = -1.0 * (group1->total_force() * dist_v.unit()); |
} else { | } else { |
group2.read_system_forces(); | group2->read_total_forces(); |
ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit()); | ft.real_value = 0.5 * ((group2->total_force() - group1->total_force()) * dist_v.unit()); |
} | } |
} | } |
| |
| |
void colvar::distance::calc_Jacobian_derivative() | void colvar::distance::calc_Jacobian_derivative() |
{ | { |
jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0; | jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0; |
} | } |
| |
| |
void colvar::distance::apply_force(colvarvalue const &force) | void colvar::distance::apply_force(colvarvalue const &force) |
{ | { |
if (!group1.noforce) | if (!group1->noforce) |
group1.apply_colvar_force(force.real_value); | group1->apply_colvar_force(force.real_value); |
| |
if (!group2.noforce) | if (!group2->noforce) |
group2.apply_colvar_force(force.real_value); | group2->apply_colvar_force(force.real_value); |
} | } |
| |
| |
| |
: distance(conf) | : distance(conf) |
{ | { |
function_type = "distance_vec"; | function_type = "distance_vec"; |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_3vector); | x.type(colvarvalue::type_3vector); |
} | } |
| |
| |
colvar::distance_vec::distance_vec() | colvar::distance_vec::distance_vec() |
: distance() | : distance() |
{ | { |
function_type = "distance_vec"; | function_type = "distance_vec"; |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_3vector); | x.type(colvarvalue::type_3vector); |
} | } |
| |
| |
void colvar::distance_vec::calc_value() | void colvar::distance_vec::calc_value() |
{ | { |
if (b_no_PBC) { | if (b_no_PBC) { |
x.rvector_value = group2.center_of_mass() - group1.center_of_mass(); | x.rvector_value = group2->center_of_mass() - group1->center_of_mass(); |
} else { | } else { |
x.rvector_value = cvm::position_distance(group1.center_of_mass(), | x.rvector_value = cvm::position_distance(group1->center_of_mass(), |
group2.center_of_mass()); | group2->center_of_mass()); |
} | } |
} | } |
| |
| |
void colvar::distance_vec::calc_gradients() | void colvar::distance_vec::calc_gradients() |
{ | { |
// gradients are not stored: a 3x3 matrix for each atom would be | // gradients are not stored: a 3x3 matrix for each atom would be |
// needed to store just the identity matrix | // needed to store just the identity matrix |
} | } |
| |
| |
void colvar::distance_vec::apply_force(colvarvalue const &force) | void colvar::distance_vec::apply_force(colvarvalue const &force) |
{ | { |
if (!group1.noforce) | if (!group1->noforce) |
group1.apply_force(-1.0 * force.rvector_value); | group1->apply_force(-1.0 * force.rvector_value); |
| |
if (!group2.noforce) | if (!group2->noforce) |
group2.apply_force( force.rvector_value); | group2->apply_force( force.rvector_value); |
} | } |
| |
| |
| |
: cvc(conf) | : cvc(conf) |
{ | { |
function_type = "distance_z"; | function_type = "distance_z"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
| |
// TODO detect PBC from MD engine (in simple cases) | // TODO detect PBC from MD engine (in simple cases) |
| |
return; | return; |
} | } |
| |
parse_group(conf, "main", main); | main = parse_group(conf, "main"); |
parse_group(conf, "ref", ref1); | ref1 = parse_group(conf, "ref"); |
atom_groups.push_back(&main); | |
atom_groups.push_back(&ref1); | |
// this group is optional | // this group is optional |
parse_group(conf, "ref2", ref2, true); | ref2 = parse_group(conf, "ref2", true); |
| |
if (ref2.size()) { | if (ref2 && ref2->size()) { |
atom_groups.push_back(&ref2); | |
cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\""); | cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\""); |
fixed_axis = false; | fixed_axis = false; |
if (key_lookup(conf, "axis")) | if (key_lookup(conf, "axis")) |
| |
if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { | if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { |
cvm::log("Computing distance using absolute positions (not minimal-image)"); | cvm::log("Computing distance using absolute positions (not minimal-image)"); |
} | } |
if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { | |
cvm::log("Computing system force on group \"main\" only"); | init_total_force_params(conf); |
} | |
} | } |
| |
colvar::distance_z::distance_z() | colvar::distance_z::distance_z() |
{ | { |
function_type = "distance_z"; | function_type = "distance_z"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
{ | { |
if (fixed_axis) { | if (fixed_axis) { |
if (b_no_PBC) { | if (b_no_PBC) { |
dist_v = main.center_of_mass() - ref1.center_of_mass(); | dist_v = main->center_of_mass() - ref1->center_of_mass(); |
} else { | } else { |
dist_v = cvm::position_distance(ref1.center_of_mass(), | dist_v = cvm::position_distance(ref1->center_of_mass(), |
main.center_of_mass()); | main->center_of_mass()); |
} | } |
} else { | } else { |
| |
if (b_no_PBC) { | if (b_no_PBC) { |
dist_v = main.center_of_mass() - | dist_v = main->center_of_mass() - |
(0.5 * (ref1.center_of_mass() + ref2.center_of_mass())); | (0.5 * (ref1->center_of_mass() + ref2->center_of_mass())); |
axis = ref2.center_of_mass() - ref1.center_of_mass(); | axis = ref2->center_of_mass() - ref1->center_of_mass(); |
} else { | } else { |
dist_v = cvm::position_distance(0.5 * (ref1.center_of_mass() + | dist_v = cvm::position_distance(0.5 * (ref1->center_of_mass() + |
ref2.center_of_mass()), main.center_of_mass()); | ref2->center_of_mass()), main->center_of_mass()); |
axis = cvm::position_distance(ref1.center_of_mass(), ref2.center_of_mass()); | axis = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass()); |
} | } |
axis_norm = axis.norm(); | axis_norm = axis.norm(); |
axis = axis.unit(); | axis = axis.unit(); |
| |
| |
void colvar::distance_z::calc_gradients() | void colvar::distance_z::calc_gradients() |
{ | { |
main.set_weighted_gradient( axis ); | main->set_weighted_gradient( axis ); |
| |
if (fixed_axis) { | if (fixed_axis) { |
ref1.set_weighted_gradient(-1.0 * axis); | ref1->set_weighted_gradient(-1.0 * axis); |
} else { | } else { |
if (b_no_PBC) { | if (b_no_PBC) { |
ref1.set_weighted_gradient( 1.0 / axis_norm * (main.center_of_mass() - ref2.center_of_mass() - | ref1->set_weighted_gradient( 1.0 / axis_norm * (main->center_of_mass() - ref2->center_of_mass() - |
x.real_value * axis )); | x.real_value * axis )); |
ref2.set_weighted_gradient( 1.0 / axis_norm * (ref1.center_of_mass() - main.center_of_mass() + | ref2->set_weighted_gradient( 1.0 / axis_norm * (ref1->center_of_mass() - main->center_of_mass() + |
x.real_value * axis )); | x.real_value * axis )); |
} else { | } else { |
ref1.set_weighted_gradient( 1.0 / axis_norm * ( | ref1->set_weighted_gradient( 1.0 / axis_norm * ( |
cvm::position_distance(ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis )); | cvm::position_distance(ref2->center_of_mass(), main->center_of_mass()) - x.real_value * axis )); |
ref2.set_weighted_gradient( 1.0 / axis_norm * ( | ref2->set_weighted_gradient( 1.0 / axis_norm * ( |
cvm::position_distance(main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis )); | cvm::position_distance(main->center_of_mass(), ref1->center_of_mass()) + x.real_value * axis )); |
} | |
} | } |
| |
if (b_debug_gradients) { | |
cvm::log("Debugging gradients for group main:\n"); | |
debug_gradients(main); | |
cvm::log("Debugging gradients for group ref1:\n"); | |
debug_gradients(ref1); | |
cvm::log("Debugging gradients for group ref2:\n"); | |
debug_gradients(ref2); | |
} | } |
} | } |
| |
void colvar::distance_z::calc_force_invgrads() | void colvar::distance_z::calc_force_invgrads() |
{ | { |
main.read_system_forces(); | main->read_total_forces(); |
| |
if (fixed_axis && !b_1site_force) { | if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) { |
ref1.read_system_forces(); | ref1->read_total_forces(); |
ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis); | ft.real_value = 0.5 * ((main->total_force() - ref1->total_force()) * axis); |
} else { | } else { |
ft.real_value = main.system_force() * axis; | ft.real_value = main->total_force() * axis; |
} | } |
} | } |
| |
| |
| |
void colvar::distance_z::apply_force(colvarvalue const &force) | void colvar::distance_z::apply_force(colvarvalue const &force) |
{ | { |
if (!ref1.noforce) | if (!ref1->noforce) |
ref1.apply_colvar_force(force.real_value); | ref1->apply_colvar_force(force.real_value); |
| |
if (ref2.size() && !ref2.noforce) | if (ref2 && ref2->size() && !ref2->noforce) |
ref2.apply_colvar_force(force.real_value); | ref2->apply_colvar_force(force.real_value); |
| |
if (!main.noforce) | if (!main->noforce) |
main.apply_colvar_force(force.real_value); | main->apply_colvar_force(force.real_value); |
} | } |
| |
| |
| |
: distance_z(conf) | : distance_z(conf) |
{ | { |
function_type = "distance_xy"; | function_type = "distance_xy"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
: distance_z() | : distance_z() |
{ | { |
function_type = "distance_xy"; | function_type = "distance_xy"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
void colvar::distance_xy::calc_value() | void colvar::distance_xy::calc_value() |
{ | { |
if (b_no_PBC) { | if (b_no_PBC) { |
dist_v = main.center_of_mass() - ref1.center_of_mass(); | dist_v = main->center_of_mass() - ref1->center_of_mass(); |
} else { | } else { |
dist_v = cvm::position_distance(ref1.center_of_mass(), | dist_v = cvm::position_distance(ref1->center_of_mass(), |
main.center_of_mass()); | main->center_of_mass()); |
} | } |
if (!fixed_axis) { | if (!fixed_axis) { |
if (b_no_PBC) { | if (b_no_PBC) { |
v12 = ref2.center_of_mass() - ref1.center_of_mass(); | v12 = ref2->center_of_mass() - ref1->center_of_mass(); |
} else { | } else { |
v12 = cvm::position_distance(ref1.center_of_mass(), ref2.center_of_mass()); | v12 = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass()); |
} | } |
axis_norm = v12.norm(); | axis_norm = v12.norm(); |
axis = v12.unit(); | axis = v12.unit(); |
| |
x_inv = 1.0 / x.real_value; | x_inv = 1.0 / x.real_value; |
| |
if (fixed_axis) { | if (fixed_axis) { |
ref1.set_weighted_gradient(-1.0 * x_inv * dist_v_ortho); | ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho); |
main.set_weighted_gradient( x_inv * dist_v_ortho); | main->set_weighted_gradient( x_inv * dist_v_ortho); |
} else { | } else { |
if (b_no_PBC) { | if (b_no_PBC) { |
v13 = main.center_of_mass() - ref1.center_of_mass(); | v13 = main->center_of_mass() - ref1->center_of_mass(); |
} else { | } else { |
v13 = cvm::position_distance(ref1.center_of_mass(), main.center_of_mass()); | v13 = cvm::position_distance(ref1->center_of_mass(), main->center_of_mass()); |
} | } |
A = (dist_v * axis) / axis_norm; | A = (dist_v * axis) / axis_norm; |
| |
ref1.set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho); | ref1->set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho); |
ref2.set_weighted_gradient( -A * x_inv * dist_v_ortho); | ref2->set_weighted_gradient( -A * x_inv * dist_v_ortho); |
main.set_weighted_gradient( 1.0 * x_inv * dist_v_ortho); | main->set_weighted_gradient( 1.0 * x_inv * dist_v_ortho); |
} | } |
} | } |
| |
void colvar::distance_xy::calc_force_invgrads() | void colvar::distance_xy::calc_force_invgrads() |
{ | { |
main.read_system_forces(); | main->read_total_forces(); |
| |
if (fixed_axis && !b_1site_force) { | if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) { |
ref1.read_system_forces(); | ref1->read_total_forces(); |
ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho); | ft.real_value = 0.5 / x.real_value * ((main->total_force() - ref1->total_force()) * dist_v_ortho); |
} else { | } else { |
ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho; | ft.real_value = 1.0 / x.real_value * main->total_force() * dist_v_ortho; |
} | } |
} | } |
| |
| |
| |
void colvar::distance_xy::apply_force(colvarvalue const &force) | void colvar::distance_xy::apply_force(colvarvalue const &force) |
{ | { |
if (!ref1.noforce) | if (!ref1->noforce) |
ref1.apply_colvar_force(force.real_value); | ref1->apply_colvar_force(force.real_value); |
| |
if (ref2.size() && !ref2.noforce) | if (ref2 && ref2->size() && !ref2->noforce) |
ref2.apply_colvar_force(force.real_value); | ref2->apply_colvar_force(force.real_value); |
| |
if (!main.noforce) | if (!main->noforce) |
main.apply_colvar_force(force.real_value); | main->apply_colvar_force(force.real_value); |
} | } |
| |
| |
| |
: distance(conf) | : distance(conf) |
{ | { |
function_type = "distance_dir"; | function_type = "distance_dir"; |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_unit3vector); | x.type(colvarvalue::type_unit3vector); |
} | } |
| |
| |
: distance() | : distance() |
{ | { |
function_type = "distance_dir"; | function_type = "distance_dir"; |
| provide(f_cvc_com_based); |
x.type(colvarvalue::type_unit3vector); | x.type(colvarvalue::type_unit3vector); |
} | } |
| |
| |
void colvar::distance_dir::calc_value() | void colvar::distance_dir::calc_value() |
{ | { |
if (b_no_PBC) { | if (b_no_PBC) { |
dist_v = group2.center_of_mass() - group1.center_of_mass(); | dist_v = group2->center_of_mass() - group1->center_of_mass(); |
} else { | } else { |
dist_v = cvm::position_distance(group1.center_of_mass(), | dist_v = cvm::position_distance(group1->center_of_mass(), |
group2.center_of_mass()); | group2->center_of_mass()); |
} | } |
x.rvector_value = dist_v.unit(); | x.rvector_value = dist_v.unit(); |
} | } |
| |
// gradients are computed on the fly within apply_force() | // gradients are computed on the fly within apply_force() |
// Note: could be a problem if a future bias relies on gradient | // Note: could be a problem if a future bias relies on gradient |
// calculations... | // calculations... |
| // TODO in new deps system: remove dependency of biasing force to gradient? |
| // That way we could tell apart an explicit gradient dependency |
} | } |
| |
| |
| |
cvm::real const iprod = force.rvector_value * x.rvector_value; | cvm::real const iprod = force.rvector_value * x.rvector_value; |
cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value; | cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value; |
| |
if (!group1.noforce) | if (!group1->noforce) |
group1.apply_force(-1.0 * force_tang); | group1->apply_force(-1.0 * force_tang); |
| |
if (!group2.noforce) | if (!group2->noforce) |
group2.apply_force( force_tang); | group2->apply_force( force_tang); |
} | } |
| |
| |
| |
return; | return; |
} | } |
| |
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { | for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { |
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { | for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { |
if (ai1->id == ai2->id) { | if (ai1->id == ai2->id) { |
cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n"); | cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n"); |
return; | return; |
| |
} | } |
} | } |
| |
b_inverse_gradients = false; | |
b_Jacobian_derivative = false; | |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
{ | { |
function_type = "distance_inv"; | function_type = "distance_inv"; |
exponent = 6; | exponent = 6; |
b_inverse_gradients = false; | |
b_Jacobian_derivative = false; | |
b_1site_force = false; | |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
{ | { |
x.real_value = 0.0; | x.real_value = 0.0; |
if (b_no_PBC) { | if (b_no_PBC) { |
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { | for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { |
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { | for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { |
cvm::rvector const dv = ai2->pos - ai1->pos; | cvm::rvector const dv = ai2->pos - ai1->pos; |
cvm::real const d2 = dv.norm2(); | cvm::real const d2 = dv.norm2(); |
cvm::real dinv = 1.0; | cvm::real dinv = 1.0; |
| |
} | } |
} | } |
} else { | } else { |
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { | for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { |
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { | for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { |
cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos); | cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos); |
cvm::real const d2 = dv.norm2(); | cvm::real const d2 = dv.norm2(); |
cvm::real dinv = 1.0; | cvm::real dinv = 1.0; |
| |
} | } |
} | } |
| |
x.real_value *= 1.0 / cvm::real(group1.size() * group2.size()); | x.real_value *= 1.0 / cvm::real(group1->size() * group2->size()); |
x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent))); | x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent))); |
} | } |
| |
void colvar::distance_inv::calc_gradients() | void colvar::distance_inv::calc_gradients() |
{ | { |
cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1.size() * group2.size()); | cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1->size() * group2->size()); |
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { | for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { |
ai1->grad *= dxdsum; | ai1->grad *= dxdsum; |
} | } |
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { | for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { |
ai2->grad *= dxdsum; | ai2->grad *= dxdsum; |
} | } |
} | } |
| |
void colvar::distance_inv::apply_force(colvarvalue const &force) | void colvar::distance_inv::apply_force(colvarvalue const &force) |
{ | { |
if (!group1.noforce) | if (!group1->noforce) |
group1.apply_colvar_force(force.real_value); | group1->apply_colvar_force(force.real_value); |
| |
if (!group2.noforce) | if (!group2->noforce) |
group2.apply_colvar_force(force.real_value); | group2->apply_colvar_force(force.real_value); |
} | } |
| |
| |
| |
: cvc(conf) | : cvc(conf) |
{ | { |
function_type = "distance_pairs"; | function_type = "distance_pairs"; |
b_inverse_gradients = false; | |
b_Jacobian_derivative = false; | |
| |
if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { | if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { |
cvm::log("Computing distance using absolute positions (not minimal-image)"); | cvm::log("Computing distance using absolute positions (not minimal-image)"); |
} | } |
| |
parse_group(conf, "group1", group1); | group1 = parse_group(conf, "group1"); |
atom_groups.push_back(&group1); | group2 = parse_group(conf, "group2"); |
| |
parse_group(conf, "group2", group2); | |
atom_groups.push_back(&group2); | |
| |
x.type(colvarvalue::type_vector); | x.type(colvarvalue::type_vector); |
x.vector1d_value.resize(group1.size() * group2.size()); | x.vector1d_value.resize(group1->size() * group2->size()); |
} | } |
| |
| |
colvar::distance_pairs::distance_pairs() | colvar::distance_pairs::distance_pairs() |
{ | { |
function_type = "distance_pairs"; | function_type = "distance_pairs"; |
b_inverse_gradients = false; | |
b_Jacobian_derivative = false; | |
x.type(colvarvalue::type_vector); | x.type(colvarvalue::type_vector); |
} | } |
| |
| |
void colvar::distance_pairs::calc_value() | void colvar::distance_pairs::calc_value() |
{ | { |
x.vector1d_value.resize(group1.size() * group2.size()); | x.vector1d_value.resize(group1->size() * group2->size()); |
| |
if (b_no_PBC) { | if (b_no_PBC) { |
size_t i1, i2; | size_t i1, i2; |
for (i1 = 0; i1 < group1.size(); i1++) { | for (i1 = 0; i1 < group1->size(); i1++) { |
for (i2 = 0; i2 < group2.size(); i2++) { | for (i2 = 0; i2 < group2->size(); i2++) { |
cvm::rvector const dv = group2[i2].pos - group1[i1].pos; | cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos; |
cvm::real const d = dv.norm(); | cvm::real const d = dv.norm(); |
x.vector1d_value[i1*group2.size() + i2] = d; | x.vector1d_value[i1*group2->size() + i2] = d; |
group1[i1].grad = -1.0 * dv.unit(); | (*group1)[i1].grad = -1.0 * dv.unit(); |
group2[i2].grad = dv.unit(); | (*group2)[i2].grad = dv.unit(); |
} | } |
} | } |
} else { | } else { |
size_t i1, i2; | size_t i1, i2; |
for (i1 = 0; i1 < group1.size(); i1++) { | for (i1 = 0; i1 < group1->size(); i1++) { |
for (i2 = 0; i2 < group2.size(); i2++) { | for (i2 = 0; i2 < group2->size(); i2++) { |
cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos); | cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos); |
cvm::real const d = dv.norm(); | cvm::real const d = dv.norm(); |
x.vector1d_value[i1*group2.size() + i2] = d; | x.vector1d_value[i1*group2->size() + i2] = d; |
group1[i1].grad = -1.0 * dv.unit(); | (*group1)[i1].grad = -1.0 * dv.unit(); |
group2[i2].grad = dv.unit(); | (*group2)[i2].grad = dv.unit(); |
} | } |
} | } |
} | } |
| |
void colvar::distance_pairs::calc_gradients() | void colvar::distance_pairs::calc_gradients() |
{ | { |
// will be calculated on the fly in apply_force() | // will be calculated on the fly in apply_force() |
if (b_debug_gradients) { | |
cvm::log("Debugging gradients:\n"); | |
debug_gradients(group1); | |
} | |
} | } |
| |
void colvar::distance_pairs::apply_force(colvarvalue const &force) | void colvar::distance_pairs::apply_force(colvarvalue const &force) |
{ | { |
if (b_no_PBC) { | if (b_no_PBC) { |
size_t i1, i2; | size_t i1, i2; |
for (i1 = 0; i1 < group1.size(); i1++) { | for (i1 = 0; i1 < group1->size(); i1++) { |
for (i2 = 0; i2 < group2.size(); i2++) { | for (i2 = 0; i2 < group2->size(); i2++) { |
cvm::rvector const dv = group2[i2].pos - group1[i1].pos; | cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos; |
group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit()); | (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit()); |
group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit()); | (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit()); |
} | } |
} | } |
} else { | } else { |
size_t i1, i2; | size_t i1, i2; |
for (i1 = 0; i1 < group1.size(); i1++) { | for (i1 = 0; i1 < group1->size(); i1++) { |
for (i2 = 0; i2 < group2.size(); i2++) { | for (i2 = 0; i2 < group2->size(); i2++) { |
cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos); | cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos); |
group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit()); | (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit()); |
group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit()); | (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit()); |
} | } |
} | } |
} | } |
| |
: cvc(conf) | : cvc(conf) |
{ | { |
function_type = "gyration"; | function_type = "gyration"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
parse_group(conf, "atoms", atoms); | atoms = parse_group(conf, "atoms"); |
atom_groups.push_back(&atoms); | |
| |
if (atoms.b_user_defined_fit) { | if (atoms->b_user_defined_fit) { |
cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); | cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); |
} else { | } else { |
atoms.b_center = true; | atoms->b_center = true; |
atoms.ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0)); | atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0)); |
} | } |
| |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
| |
colvar::gyration::gyration() | colvar::gyration::gyration() |
{ | { |
function_type = "gyration"; | function_type = "gyration"; |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
void colvar::gyration::calc_value() | void colvar::gyration::calc_value() |
{ | { |
x.real_value = 0.0; | x.real_value = 0.0; |
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { | for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { |
x.real_value += (ai->pos).norm2(); | x.real_value += (ai->pos).norm2(); |
} | } |
x.real_value = std::sqrt(x.real_value / cvm::real(atoms.size())); | x.real_value = std::sqrt(x.real_value / cvm::real(atoms->size())); |
} | } |
| |
| |
void colvar::gyration::calc_gradients() | void colvar::gyration::calc_gradients() |
{ | { |
cvm::real const drdx = 1.0/(cvm::real(atoms.size()) * x.real_value); | cvm::real const drdx = 1.0/(cvm::real(atoms->size()) * x.real_value); |
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { | for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { |
ai->grad = drdx * ai->pos; | ai->grad = drdx * ai->pos; |
} | } |
| |
if (b_debug_gradients) { | |
cvm::log("Debugging gradients:\n"); | |
debug_gradients(atoms); | |
} | |
} | } |
| |
| |
void colvar::gyration::calc_force_invgrads() | void colvar::gyration::calc_force_invgrads() |
{ | { |
atoms.read_system_forces(); | atoms->read_total_forces(); |
| |
cvm::real const dxdr = 1.0/x.real_value; | cvm::real const dxdr = 1.0/x.real_value; |
ft.real_value = 0.0; | ft.real_value = 0.0; |
| |
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { | for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { |
ft.real_value += dxdr * ai->pos * ai->system_force; | ft.real_value += dxdr * ai->pos * ai->total_force; |
} | } |
} | } |
| |
| |
void colvar::gyration::calc_Jacobian_derivative() | void colvar::gyration::calc_Jacobian_derivative() |
{ | { |
jd = x.real_value ? (3.0 * cvm::real(atoms.size()) - 4.0) / x.real_value : 0.0; | jd = x.real_value ? (3.0 * cvm::real(atoms->size()) - 4.0) / x.real_value : 0.0; |
} | } |
| |
| |
void colvar::gyration::apply_force(colvarvalue const &force) | void colvar::gyration::apply_force(colvarvalue const &force) |
{ | { |
if (!atoms.noforce) | if (!atoms->noforce) |
atoms.apply_colvar_force(force.real_value); | atoms->apply_colvar_force(force.real_value); |
} | } |
| |
| |
| |
: gyration(conf) | : gyration(conf) |
{ | { |
function_type = "inertia"; | function_type = "inertia"; |
b_inverse_gradients = false; | |
b_Jacobian_derivative = false; | |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
} | } |
| |
| |
void colvar::inertia::calc_value() | void colvar::inertia::calc_value() |
{ | { |
x.real_value = 0.0; | x.real_value = 0.0; |
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { | for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { |
x.real_value += (ai->pos).norm2(); | x.real_value += (ai->pos).norm2(); |
} | } |
} | } |
| |
| |
void colvar::inertia::calc_gradients() | void colvar::inertia::calc_gradients() |
{ | { |
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { | for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { |
ai->grad = 2.0 * ai->pos; | ai->grad = 2.0 * ai->pos; |
} | } |
| |
if (b_debug_gradients) { | |
cvm::log("Debugging gradients:\n"); | |
debug_gradients(atoms); | |
} | |
} | } |
| |
| |
void colvar::inertia::apply_force(colvarvalue const &force) | void colvar::inertia::apply_force(colvarvalue const &force) |
{ | { |
if (!atoms.noforce) | if (!atoms->noforce) |
atoms.apply_colvar_force(force.real_value); | atoms->apply_colvar_force(force.real_value); |
} | } |
| |
| |
| |
void colvar::inertia_z::calc_value() | void colvar::inertia_z::calc_value() |
{ | { |
x.real_value = 0.0; | x.real_value = 0.0; |
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { | for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { |
cvm::real const iprod = ai->pos * axis; | cvm::real const iprod = ai->pos * axis; |
x.real_value += iprod * iprod; | x.real_value += iprod * iprod; |
} | } |
| |
| |
void colvar::inertia_z::calc_gradients() | void colvar::inertia_z::calc_gradients() |
{ | { |
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { | for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { |
ai->grad = 2.0 * (ai->pos * axis) * axis; | ai->grad = 2.0 * (ai->pos * axis) * axis; |
} | } |
| |
if (b_debug_gradients) { | |
cvm::log("Debugging gradients:\n"); | |
debug_gradients(atoms); | |
} | |
} | } |
| |
| |
void colvar::inertia_z::apply_force(colvarvalue const &force) | void colvar::inertia_z::apply_force(colvarvalue const &force) |
{ | { |
if (!atoms.noforce) | if (!atoms->noforce) |
atoms.apply_colvar_force(force.real_value); | atoms->apply_colvar_force(force.real_value); |
} | } |
| |
| |
| |
colvar::rmsd::rmsd(std::string const &conf) | colvar::rmsd::rmsd(std::string const &conf) |
: cvc(conf) | : cvc(conf) |
{ | { |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | |
function_type = "rmsd"; | function_type = "rmsd"; |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
| |
parse_group(conf, "atoms", atoms); | atoms = parse_group(conf, "atoms"); |
atom_groups.push_back(&atoms); | |
| |
if (atoms.b_dummy) { | if (!atoms || atoms->size() == 0) { |
cvm::error("Error: \"atoms\" cannot be a dummy atom."); | cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD."); |
return; | return; |
} | } |
| |
if (atoms.ref_pos_group != NULL && b_Jacobian_derivative) { | bool b_Jacobian_derivative = true; |
| if (atoms->fitting_group != NULL && b_Jacobian_derivative) { |
cvm::log("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: " | cvm::log("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: " |
"Jacobian derivatives of the RMSD will not be calculated.\n"); | "Jacobian derivatives of the RMSD will not be calculated.\n"); |
b_Jacobian_derivative = false; | b_Jacobian_derivative = false; |
} | } |
| if (b_Jacobian_derivative) provide(f_cvc_Jacobian); |
| |
// the following is a simplified version of the corresponding atom group options; | // the following is a simplified version of the corresponding atom group options; |
// we need this because the reference coordinates defined inside the atom group | // we need this because the reference coordinates defined inside the atom group |
// may be used only for fitting, and even more so if ref_pos_group is used | // may be used only for fitting, and even more so if fitting_group is used |
if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) { | if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) { |
cvm::log("Using reference positions from configuration file to calculate the variable.\n"); | cvm::log("Using reference positions from configuration file to calculate the variable.\n"); |
if (ref_pos.size() != atoms.size()) { | if (ref_pos.size() != atoms->size()) { |
cvm::error("Error: the number of reference positions provided ("+ | cvm::error("Error: the number of reference positions provided ("+ |
cvm::to_str(ref_pos.size())+ | cvm::to_str(ref_pos.size())+ |
") does not match the number of atoms of group \"atoms\" ("+ | ") does not match the number of atoms of group \"atoms\" ("+ |
cvm::to_str(atoms.size())+").\n"); | cvm::to_str(atoms->size())+").\n"); |
return; | return; |
} | } |
} | } |
| |
} | } |
} else { | } else { |
// if not, rely on existing atom indices for the group | // if not, rely on existing atom indices for the group |
atoms.create_sorted_ids(); | atoms->create_sorted_ids(); |
| ref_pos.resize(atoms->size()); |
} | } |
| |
ref_pos.resize(atoms.size()); | cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms->sorted_ids, |
cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms.sorted_ids, | |
ref_pos_col, ref_pos_col_value); | ref_pos_col, ref_pos_col_value); |
} | } |
} | } |
| |
if (atoms.b_user_defined_fit) { | if (ref_pos.size() != atoms->size()) { |
| cvm::error("Error: found " + cvm::to_str(ref_pos.size()) + |
| " reference positions; expected " + cvm::to_str(atoms->size())); |
| return; |
| } |
| |
| if (atoms->b_user_defined_fit) { |
cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); | cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); |
} else { | } else { |
// Default: fit everything | // Default: fit everything |
cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: " | cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: " |
"if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n"); | "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n"); |
atoms.b_center = true; | atoms->b_center = true; |
atoms.b_rotate = true; | atoms->b_rotate = true; |
// default case: reference positions for calculating the rmsd are also those used | // default case: reference positions for calculating the rmsd are also those used |
// for fitting | // for fitting |
atoms.ref_pos = ref_pos; | atoms->ref_pos = ref_pos; |
atoms.center_ref_pos(); | atoms->center_ref_pos(); |
| |
cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation " | cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation " |
"will not be computed as they cancel out in the gradients."); | "will not be computed as they cancel out in the gradients."); |
atoms.b_fit_gradients = false; | atoms->b_fit_gradients = false; |
} | |
| |
if (atoms.b_rotate) { | |
// TODO: finer-grained control of this would require exposing a | |
// "request_Jacobian_derivative()" method to the colvar, and the same | |
// from the colvar to biases | |
// TODO: this should not be enabled here anyway, as it is not specific of the | |
// component - instead it should be decided in a generic way by the atom group | |
| |
// request the calculation of the derivatives of the rotation defined by the atom group | // request the calculation of the derivatives of the rotation defined by the atom group |
atoms.rot.request_group1_gradients(atoms.size()); | atoms->rot.request_group1_gradients(atoms->size()); |
// request derivatives of optimal rotation wrt reference coordinates for Jacobian: | // request derivatives of optimal rotation wrt reference coordinates for Jacobian: |
atoms.rot.request_group2_gradients(atoms.size()); | // this is only required for ABF, but we do both groups here for better caching |
| atoms->rot.request_group2_gradients(atoms->size()); |
} | } |
} | } |
| |
| |
// rotational-translational fit is handled by the atom group | // rotational-translational fit is handled by the atom group |
| |
x.real_value = 0.0; | x.real_value = 0.0; |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2(); | x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2(); |
} | } |
x.real_value /= cvm::real(atoms.size()); // MSD | x.real_value /= cvm::real(atoms->size()); // MSD |
x.real_value = std::sqrt(x.real_value); | x.real_value = std::sqrt(x.real_value); |
} | } |
| |
| |
void colvar::rmsd::calc_gradients() | void colvar::rmsd::calc_gradients() |
{ | { |
cvm::real const drmsddx2 = (x.real_value > 0.0) ? | cvm::real const drmsddx2 = (x.real_value > 0.0) ? |
0.5 / (x.real_value * cvm::real(atoms.size())) : | 0.5 / (x.real_value * cvm::real(atoms->size())) : |
0.0; | 0.0; |
| |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - ref_pos[ia])); | (*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[ia])); |
} | |
| |
if (b_debug_gradients) { | |
cvm::log("Debugging gradients:\n"); | |
debug_gradients(atoms); | |
} | } |
} | } |
| |
| |
void colvar::rmsd::apply_force(colvarvalue const &force) | void colvar::rmsd::apply_force(colvarvalue const &force) |
{ | { |
if (!atoms.noforce) | if (!atoms->noforce) |
atoms.apply_colvar_force(force.real_value); | atoms->apply_colvar_force(force.real_value); |
} | } |
| |
| |
void colvar::rmsd::calc_force_invgrads() | void colvar::rmsd::calc_force_invgrads() |
{ | { |
atoms.read_system_forces(); | atoms->read_total_forces(); |
ft.real_value = 0.0; | ft.real_value = 0.0; |
| |
// Note: gradient square norm is 1/N_atoms | // Note: gradient square norm is 1/N_atoms |
| |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
ft.real_value += atoms[ia].grad * atoms[ia].system_force; | ft.real_value += (*atoms)[ia].grad * (*atoms)[ia].total_force; |
} | } |
ft.real_value *= atoms.size(); | ft.real_value *= atoms->size(); |
} | } |
| |
| |
| |
// divergence of the rotated coordinates (including only derivatives of the rotation matrix) | // divergence of the rotated coordinates (including only derivatives of the rotation matrix) |
cvm::real divergence = 0.0; | cvm::real divergence = 0.0; |
| |
if (atoms.b_rotate) { | if (atoms->b_rotate) { |
| |
// gradient of the rotation matrix | // gradient of the rotation matrix |
cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3); | cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3); |
// gradients of products of 2 quaternion components | // gradients of products of 2 quaternion components |
cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; | cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
| |
// Gradient of optimal quaternion wrt current Cartesian position | // Gradient of optimal quaternion wrt current Cartesian position |
cvm::vector1d<cvm::rvector> &dq = atoms.rot.dQ0_1[ia]; | cvm::vector1d<cvm::rvector> &dq = atoms->rot.dQ0_1[ia]; |
| |
g11 = 2.0 * (atoms.rot.q)[1]*dq[1]; | g11 = 2.0 * (atoms->rot.q)[1]*dq[1]; |
g22 = 2.0 * (atoms.rot.q)[2]*dq[2]; | g22 = 2.0 * (atoms->rot.q)[2]*dq[2]; |
g33 = 2.0 * (atoms.rot.q)[3]*dq[3]; | g33 = 2.0 * (atoms->rot.q)[3]*dq[3]; |
g01 = (atoms.rot.q)[0]*dq[1] + (atoms.rot.q)[1]*dq[0]; | g01 = (atoms->rot.q)[0]*dq[1] + (atoms->rot.q)[1]*dq[0]; |
g02 = (atoms.rot.q)[0]*dq[2] + (atoms.rot.q)[2]*dq[0]; | g02 = (atoms->rot.q)[0]*dq[2] + (atoms->rot.q)[2]*dq[0]; |
g03 = (atoms.rot.q)[0]*dq[3] + (atoms.rot.q)[3]*dq[0]; | g03 = (atoms->rot.q)[0]*dq[3] + (atoms->rot.q)[3]*dq[0]; |
g12 = (atoms.rot.q)[1]*dq[2] + (atoms.rot.q)[2]*dq[1]; | g12 = (atoms->rot.q)[1]*dq[2] + (atoms->rot.q)[2]*dq[1]; |
g13 = (atoms.rot.q)[1]*dq[3] + (atoms.rot.q)[3]*dq[1]; | g13 = (atoms->rot.q)[1]*dq[3] + (atoms->rot.q)[3]*dq[1]; |
g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2]; | g23 = (atoms->rot.q)[2]*dq[3] + (atoms->rot.q)[3]*dq[2]; |
| |
// Gradient of the rotation matrix wrt current Cartesian position | // Gradient of the rotation matrix wrt current Cartesian position |
grad_rot_mat[0][0] = -2.0 * (g22 + g33); | grad_rot_mat[0][0] = -2.0 * (g22 + g33); |
| |
} | } |
} | } |
} | } |
| jd.real_value = x.real_value > 0.0 ? (3.0 * atoms->size() - 4.0 - divergence) / x.real_value : 0.0; |
jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0; | |
} | } |
| |
| |
| |
colvar::eigenvector::eigenvector(std::string const &conf) | colvar::eigenvector::eigenvector(std::string const &conf) |
: cvc(conf) | : cvc(conf) |
{ | { |
b_inverse_gradients = true; | provide(f_cvc_inv_gradient); |
b_Jacobian_derivative = true; | provide(f_cvc_Jacobian); |
function_type = "eigenvector"; | function_type = "eigenvector"; |
x.type(colvarvalue::type_scalar); | x.type(colvarvalue::type_scalar); |
| |
parse_group(conf, "atoms", atoms); | atoms = parse_group(conf, "atoms"); |
atom_groups.push_back(&atoms); | |
| |
| |
{ | { |
bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos); | bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos); |
| |
if (b_inline) { | if (b_inline) { |
cvm::log("Using reference positions from input file.\n"); | cvm::log("Using reference positions from input file.\n"); |
if (ref_pos.size() != atoms.size()) { | if (ref_pos.size() != atoms->size()) { |
cvm::error("Error: reference positions do not " | cvm::error("Error: reference positions do not " |
"match the number of requested atoms.\n"); | "match the number of requested atoms.\n"); |
return; | return; |
| |
} | } |
} else { | } else { |
// if not, use atom indices | // if not, use atom indices |
atoms.create_sorted_ids(); | atoms->create_sorted_ids(); |
} | } |
| |
ref_pos.resize(atoms.size()); | ref_pos.resize(atoms->size()); |
cvm::load_coords(file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value); | cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value); |
} | } |
} | } |
| |
| if (ref_pos.size() == 0) { |
| cvm::error("Error: reference positions were not provided.\n", INPUT_ERROR); |
| return; |
| } |
| |
| if (ref_pos.size() != atoms->size()) { |
| cvm::error("Error: reference positions do not " |
| "match the number of requested atoms.\n", INPUT_ERROR); |
| return; |
| } |
| |
// save for later the geometric center of the provided positions (may not be the origin) | // save for later the geometric center of the provided positions (may not be the origin) |
cvm::rvector ref_pos_center(0.0, 0.0, 0.0); | cvm::rvector ref_pos_center(0.0, 0.0, 0.0); |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
ref_pos_center += ref_pos[i]; | ref_pos_center += ref_pos[i]; |
} | } |
ref_pos_center *= 1.0 / atoms.size(); | ref_pos_center *= 1.0 / atoms->size(); |
| |
if (atoms.b_user_defined_fit) { | if (atoms->b_user_defined_fit) { |
cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n"); | cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n"); |
} else { | } else { |
// default: fit everything | // default: fit everything |
cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: " | cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: " |
"if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n"); | "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n"); |
atoms.b_center = true; | atoms->b_center = true; |
atoms.b_rotate = true; | atoms->b_rotate = true; |
atoms.ref_pos = ref_pos; | atoms->ref_pos = ref_pos; |
atoms.center_ref_pos(); | atoms->center_ref_pos(); |
| atoms->b_fit_gradients = false; // cancel out if group is fitted on itself |
| // and cvc is translationally invariant |
| |
// request the calculation of the derivatives of the rotation defined by the atom group | // request the calculation of the derivatives of the rotation defined by the atom group |
atoms.rot.request_group1_gradients(atoms.size()); | atoms->rot.request_group1_gradients(atoms->size()); |
// request derivatives of optimal rotation wrt reference coordinates for Jacobian: | // request derivatives of optimal rotation wrt reference coordinates for Jacobian: |
// this is only required for ABF, but we do both groups here for better caching | // this is only required for ABF, but we do both groups here for better caching |
atoms.rot.request_group2_gradients(atoms.size()); | atoms->rot.request_group2_gradients(atoms->size()); |
} | } |
| |
{ | { |
| |
// now load the eigenvector | // now load the eigenvector |
if (b_inline) { | if (b_inline) { |
cvm::log("Using vector components from input file.\n"); | cvm::log("Using vector components from input file.\n"); |
if (eigenvec.size() != atoms.size()) { | if (eigenvec.size() != atoms->size()) { |
cvm::fatal_error("Error: vector components do not " | cvm::fatal_error("Error: vector components do not " |
"match the number of requested atoms.\n"); | "match the number of requested atoms->\n"); |
} | } |
} | } |
| |
| |
} | } |
} else { | } else { |
// if not, use atom indices | // if not, use atom indices |
atoms.create_sorted_ids(); | atoms->create_sorted_ids(); |
} | } |
| |
eigenvec.resize(atoms.size()); | eigenvec.resize(atoms->size()); |
cvm::load_coords(file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value); | cvm::load_coords(file_name.c_str(), eigenvec, atoms->sorted_ids, file_col, file_col_value); |
} | } |
} | } |
| |
| |
} | } |
| |
cvm::atom_pos eig_center(0.0, 0.0, 0.0); | cvm::atom_pos eig_center(0.0, 0.0, 0.0); |
for (size_t eil = 0; eil < atoms.size(); eil++) { | for (size_t eil = 0; eil < atoms->size(); eil++) { |
eig_center += eigenvec[eil]; | eig_center += eigenvec[eil]; |
} | } |
eig_center *= 1.0 / atoms.size(); | eig_center *= 1.0 / atoms->size(); |
cvm::log("Geometric center of the provided vector: "+cvm::to_str(eig_center)+"\n"); | cvm::log("Geometric center of the provided vector: "+cvm::to_str(eig_center)+"\n"); |
| |
bool b_difference_vector = false; | bool b_difference_vector = false; |
| |
| |
if (b_difference_vector) { | if (b_difference_vector) { |
| |
if (atoms.b_center) { | if (atoms->b_center) { |
// both sets should be centered on the origin for fitting | // both sets should be centered on the origin for fitting |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
eigenvec[i] -= eig_center; | eigenvec[i] -= eig_center; |
ref_pos[i] -= ref_pos_center; | ref_pos[i] -= ref_pos_center; |
} | } |
} | } |
if (atoms.b_rotate) { | if (atoms->b_rotate) { |
atoms.rot.calc_optimal_rotation(eigenvec, ref_pos); | atoms->rot.calc_optimal_rotation(eigenvec, ref_pos); |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
eigenvec[i] = atoms.rot.rotate(eigenvec[i]); | eigenvec[i] = atoms->rot.rotate(eigenvec[i]); |
} | } |
} | } |
cvm::log("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n"); | cvm::log("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n"); |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
eigenvec[i] -= ref_pos[i]; | eigenvec[i] -= ref_pos[i]; |
} | } |
if (atoms.b_center) { | if (atoms->b_center) { |
// bring back the ref positions to where they were | // bring back the ref positions to where they were |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
ref_pos[i] += ref_pos_center; | ref_pos[i] += ref_pos_center; |
} | } |
} | } |
| |
} else { | } else { |
cvm::log("Centering the provided vector to zero.\n"); | cvm::log("Centering the provided vector to zero.\n"); |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
eigenvec[i] -= eig_center; | eigenvec[i] -= eig_center; |
} | } |
} | } |
| |
| |
// for inverse gradients | // for inverse gradients |
eigenvec_invnorm2 = 0.0; | eigenvec_invnorm2 = 0.0; |
for (size_t ein = 0; ein < atoms.size(); ein++) { | for (size_t ein = 0; ein < atoms->size(); ein++) { |
eigenvec_invnorm2 += eigenvec[ein].norm2(); | eigenvec_invnorm2 += eigenvec[ein].norm2(); |
} | } |
eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2; | eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2; |
| |
if (b_difference_vector) { | if (b_difference_vector) { |
cvm::log("\"differenceVector\" is on: normalizing the vector.\n"); | cvm::log("\"differenceVector\" is on: normalizing the vector.\n"); |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
eigenvec[i] *= eigenvec_invnorm2; | eigenvec[i] *= eigenvec_invnorm2; |
} | } |
} else { | } else { |
| |
void colvar::eigenvector::calc_value() | void colvar::eigenvector::calc_value() |
{ | { |
x.real_value = 0.0; | x.real_value = 0.0; |
for (size_t i = 0; i < atoms.size(); i++) { | for (size_t i = 0; i < atoms->size(); i++) { |
x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i]; | x.real_value += ((*atoms)[i].pos - ref_pos[i]) * eigenvec[i]; |
} | } |
} | } |
| |
| |
void colvar::eigenvector::calc_gradients() | void colvar::eigenvector::calc_gradients() |
{ | { |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
atoms[ia].grad = eigenvec[ia]; | (*atoms)[ia].grad = eigenvec[ia]; |
} | |
| |
if (b_debug_gradients) { | |
cvm::log("Debugging gradients:\n"); | |
debug_gradients(atoms); | |
} | } |
} | } |
| |
| |
void colvar::eigenvector::apply_force(colvarvalue const &force) | void colvar::eigenvector::apply_force(colvarvalue const &force) |
{ | { |
if (!atoms.noforce) | if (!atoms->noforce) |
atoms.apply_colvar_force(force.real_value); | atoms->apply_colvar_force(force.real_value); |
} | } |
| |
| |
void colvar::eigenvector::calc_force_invgrads() | void colvar::eigenvector::calc_force_invgrads() |
{ | { |
atoms.read_system_forces(); | atoms->read_total_forces(); |
ft.real_value = 0.0; | ft.real_value = 0.0; |
| |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
ft.real_value += eigenvec_invnorm2 * atoms[ia].grad * | ft.real_value += eigenvec_invnorm2 * (*atoms)[ia].grad * |
atoms[ia].system_force; | (*atoms)[ia].total_force; |
} | } |
} | } |
| |
| |
{ | { |
// gradient of the rotation matrix | // gradient of the rotation matrix |
cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3); | cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3); |
cvm::quaternion &quat0 = atoms.rot.q; | cvm::quaternion &quat0 = atoms->rot.q; |
| |
// gradients of products of 2 quaternion components | // gradients of products of 2 quaternion components |
cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; | cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; |
| |
cvm::real sum = 0.0; | cvm::real sum = 0.0; |
| |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
| |
// Gradient of optimal quaternion wrt current Cartesian position | // Gradient of optimal quaternion wrt current Cartesian position |
// trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t | // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t |
// we can just transpose the derivatives of the direct matrix | // we can just transpose the derivatives of the direct matrix |
cvm::vector1d<cvm::rvector> &dq_1 = atoms.rot.dQ0_1[ia]; | cvm::vector1d<cvm::rvector> &dq_1 = atoms->rot.dQ0_1[ia]; |
| |
g11 = 2.0 * quat0[1]*dq_1[1]; | g11 = 2.0 * quat0[1]*dq_1[1]; |
g22 = 2.0 * quat0[2]*dq_1[2]; | g22 = 2.0 * quat0[2]*dq_1[2]; |
| |
colvar::cartesian::cartesian(std::string const &conf) | colvar::cartesian::cartesian(std::string const &conf) |
: cvc(conf) | : cvc(conf) |
{ | { |
b_inverse_gradients = false; | |
b_Jacobian_derivative = false; | |
function_type = "cartesian"; | function_type = "cartesian"; |
| |
parse_group(conf, "atoms", atoms); | atoms = parse_group(conf, "atoms"); |
atom_groups.push_back(&atoms); | |
| |
bool use_x, use_y, use_z; | bool use_x, use_y, use_z; |
get_keyval(conf, "useX", use_x, true); | get_keyval(conf, "useX", use_x, true); |
| |
| |
if (axes.size() == 0) { | if (axes.size() == 0) { |
cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n"); | cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n"); |
| return; |
} | } |
| |
x.type(colvarvalue::type_vector); | x.type(colvarvalue::type_vector); |
x.vector1d_value.resize(atoms.size() * axes.size()); | x.vector1d_value.resize(atoms->size() * axes.size()); |
} | } |
| |
| |
| |
{ | { |
size_t const dim = axes.size(); | size_t const dim = axes.size(); |
size_t ia, j; | size_t ia, j; |
for (ia = 0; ia < atoms.size(); ia++) { | for (ia = 0; ia < atoms->size(); ia++) { |
for (j = 0; j < dim; j++) { | for (j = 0; j < dim; j++) { |
x.vector1d_value[dim*ia + j] = atoms[ia].pos[axes[j]]; | x.vector1d_value[dim*ia + j] = (*atoms)[ia].pos[axes[j]]; |
} | } |
} | } |
} | } |
| |
{ | { |
size_t const dim = axes.size(); | size_t const dim = axes.size(); |
size_t ia, j; | size_t ia, j; |
if (!atoms.noforce) { | if (!atoms->noforce) { |
cvm::rvector f; | cvm::rvector f; |
for (ia = 0; ia < atoms.size(); ia++) { | for (ia = 0; ia < atoms->size(); ia++) { |
for (j = 0; j < dim; j++) { | for (j = 0; j < dim; j++) { |
f[axes[j]] = force.vector1d_value[dim*ia + j]; | f[axes[j]] = force.vector1d_value[dim*ia + j]; |
} | } |
atoms[ia].apply_force(f); | (*atoms)[ia].apply_force(f); |
} | } |
} | } |
} | } |