version 1.15 | version 1.16 |
---|
| |
/// -*- c++ -*- | // -*- c++ -*- |
| |
#include <cmath> | #include <cmath> |
| |
| |
: cvc(conf) | : cvc(conf) |
{ | { |
function_type = "orientation"; | function_type = "orientation"; |
parse_group(conf, "atoms", atoms); | atoms = parse_group(conf, "atoms"); |
atom_groups.push_back(&atoms); | |
x.type(colvarvalue::type_quaternion); | x.type(colvarvalue::type_quaternion); |
| |
ref_pos.reserve(atoms.size()); | ref_pos.reserve(atoms->size()); |
| |
if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) { | if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) { |
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::fatal_error("Error: reference positions do not " | cvm::fatal_error("Error: reference positions do not " |
"match the number of requested atoms.\n"); | "match the number of requested atoms.\n"); |
} | } |
| |
"if provided, must be non-zero.\n"); | "if provided, must be non-zero.\n"); |
} 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); |
} | } |
} | } |
| |
| |
cvm::log("Centering the reference coordinates: it is " | cvm::log("Centering the reference coordinates: it is " |
"assumed that each atom is the closest " | "assumed that each atom is the closest " |
"periodic image to the center of geometry.\n"); | "periodic image to the center of geometry.\n"); |
cvm::rvector cog(0.0, 0.0, 0.0); | cvm::rvector ref_cog(0.0, 0.0, 0.0); |
size_t i; | size_t i; |
for (i = 0; i < ref_pos.size(); i++) { | for (i = 0; i < ref_pos.size(); i++) { |
cog += ref_pos[i]; | ref_cog += ref_pos[i]; |
} | } |
cog /= cvm::real(ref_pos.size()); | ref_cog /= cvm::real(ref_pos.size()); |
for (i = 0; i < ref_pos.size(); i++) { | for (i = 0; i < ref_pos.size(); i++) { |
ref_pos[i] -= cog; | ref_pos[i] -= ref_cog; |
} | } |
| |
get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0)); | get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0)); |
| |
// initialize rot member data | // initialize rot member data |
if (!atoms.noforce) { | if (!atoms->noforce) { |
rot.request_group2_gradients(atoms.size()); | rot.request_group2_gradients(atoms->size()); |
} | } |
| |
} | } |
| |
| |
void colvar::orientation::calc_value() | void colvar::orientation::calc_value() |
{ | { |
// atoms.reset_atoms_data(); | rot.b_debug_gradients = is_enabled(f_cvc_debug_gradient); |
// atoms.read_positions(); | atoms_cog = atoms->center_of_geometry(); |
| |
atoms_cog = atoms.center_of_geometry(); | rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); |
| |
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); | |
| |
if ((rot.q).inner(ref_quat) >= 0.0) { | if ((rot.q).inner(ref_quat) >= 0.0) { |
x.quaternion_value = rot.q; | x.quaternion_value = rot.q; |
| |
{ | { |
cvm::quaternion const &FQ = force.quaternion_value; | cvm::quaternion const &FQ = force.quaternion_value; |
| |
if (!atoms.noforce) { | if (!atoms->noforce) { |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
for (size_t i = 0; i < 4; i++) { | for (size_t i = 0; i < 4; i++) { |
atoms[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]); | (*atoms)[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]); |
} | } |
} | } |
} | } |
| |
| |
void colvar::orientation_angle::calc_value() | void colvar::orientation_angle::calc_value() |
{ | { |
// atoms.reset_atoms_data(); | atoms_cog = atoms->center_of_geometry(); |
// atoms.read_positions(); | |
| |
atoms_cog = atoms.center_of_geometry(); | |
| |
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); | rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); |
| |
if ((rot.q).q0 >= 0.0) { | if ((rot.q).q0 >= 0.0) { |
x.real_value = (180.0/PI) * 2.0 * std::acos((rot.q).q0); | x.real_value = (180.0/PI) * 2.0 * std::acos((rot.q).q0); |
| |
((180.0 / PI) * (-2.0) / std::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) : | ((180.0 / PI) * (-2.0) / std::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) : |
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 = (dxdq0 * (rot.dQ0_2[ia])[0]); | (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); |
} | |
if (b_debug_gradients) { | |
cvm::log("Debugging orientationAngle component gradients:\n"); | |
debug_gradients(atoms); | |
} | } |
} | } |
| |
| |
void colvar::orientation_angle::apply_force(colvarvalue const &force) | void colvar::orientation_angle::apply_force(colvarvalue const &force) |
{ | { |
cvm::real const &fw = force.real_value; | cvm::real const &fw = force.real_value; |
if (!atoms.noforce) { | if (!atoms->noforce) { |
atoms.apply_colvar_force(fw); | atoms->apply_colvar_force(fw); |
} | } |
} | } |
| |
| |
| |
void colvar::orientation_proj::calc_value() | void colvar::orientation_proj::calc_value() |
{ | { |
atoms_cog = atoms.center_of_geometry(); | atoms_cog = atoms->center_of_geometry(); |
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); | rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); |
x.real_value = 2.0 * (rot.q).q0 * (rot.q).q0 - 1.0; | x.real_value = 2.0 * (rot.q).q0 * (rot.q).q0 - 1.0; |
} | } |
| |
| |
void colvar::orientation_proj::calc_gradients() | void colvar::orientation_proj::calc_gradients() |
{ | { |
cvm::real const dxdq0 = 2.0 * 2.0 * (rot.q).q0; | cvm::real const dxdq0 = 2.0 * 2.0 * (rot.q).q0; |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); | (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); |
} | |
if (b_debug_gradients) { | |
cvm::log("Debugging orientationProj component gradients:\n"); | |
debug_gradients(atoms); | |
} | } |
} | } |
| |
| |
{ | { |
cvm::real const &fw = force.real_value; | cvm::real const &fw = force.real_value; |
| |
if (!atoms.noforce) { | if (!atoms->noforce) { |
atoms.apply_colvar_force(fw); | atoms->apply_colvar_force(fw); |
} | } |
} | } |
| |
| |
| |
void colvar::tilt::calc_value() | void colvar::tilt::calc_value() |
{ | { |
// atoms.reset_atoms_data(); | atoms_cog = atoms->center_of_geometry(); |
// atoms.read_positions(); | |
| |
atoms_cog = atoms.center_of_geometry(); | rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); |
| |
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); | |
| |
x.real_value = rot.cos_theta(axis); | x.real_value = rot.cos_theta(axis); |
} | } |
| |
{ | { |
cvm::quaternion const dxdq = rot.dcos_theta_dq(axis); | cvm::quaternion const dxdq = rot.dcos_theta_dq(axis); |
| |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
atoms[ia].grad = cvm::rvector(0.0, 0.0, 0.0); | (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0); |
for (size_t iq = 0; iq < 4; iq++) { | for (size_t iq = 0; iq < 4; iq++) { |
atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); | (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); |
} | |
} | } |
| |
if (b_debug_gradients) { | |
cvm::log("Debugging tilt component gradients:\n"); | |
debug_gradients(atoms); | |
} | } |
} | } |
| |
| |
{ | { |
cvm::real const &fw = force.real_value; | cvm::real const &fw = force.real_value; |
| |
if (!atoms.noforce) { | if (!atoms->noforce) { |
atoms.apply_colvar_force(fw); | atoms->apply_colvar_force(fw); |
} | } |
} | } |
| |
| |
| |
void colvar::spin_angle::calc_value() | void colvar::spin_angle::calc_value() |
{ | { |
// atoms.reset_atoms_data(); | atoms_cog = atoms->center_of_geometry(); |
// atoms.read_positions(); | |
| |
atoms_cog = atoms.center_of_geometry(); | |
| |
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); | rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); |
| |
x.real_value = rot.spin_angle(axis); | x.real_value = rot.spin_angle(axis); |
this->wrap(x); | this->wrap(x); |
| |
{ | { |
cvm::quaternion const dxdq = rot.dspin_angle_dq(axis); | cvm::quaternion const dxdq = rot.dspin_angle_dq(axis); |
| |
for (size_t ia = 0; ia < atoms.size(); ia++) { | for (size_t ia = 0; ia < atoms->size(); ia++) { |
atoms[ia].grad = cvm::rvector(0.0, 0.0, 0.0); | (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0); |
for (size_t iq = 0; iq < 4; iq++) { | for (size_t iq = 0; iq < 4; iq++) { |
atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); | (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); |
} | } |
} | } |
} | } |
| |
{ | { |
cvm::real const &fw = force.real_value; | cvm::real const &fw = force.real_value; |
| |
if (!atoms.noforce) { | if (!atoms->noforce) { |
atoms.apply_colvar_force(fw); | atoms->apply_colvar_force(fw); |
} | } |
} | } |