Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

colvar::eigenvector Class Reference

Colvar component: projection of 3N coordinates onto an eigenvector (colvarvalue::type_scalar type, range (-*:*)). More...

#include <colvarcomp.h>

Inheritance diagram for colvar::eigenvector:

colvar::cvc colvarparse List of all members.

Public Member Functions

 eigenvector (std::string const &conf)
 Constructor.
virtual ~eigenvector ()
virtual void calc_value ()
 Calculate the variable.
virtual void calc_gradients ()
 Calculate the atomic gradients, to be reused later in order to apply forces.
virtual void calc_force_invgrads ()
 Calculate the total force from the system using the inverse atomic gradients.
virtual void calc_Jacobian_derivative ()
 Calculate the divergence of the inverse atomic gradients.
virtual void apply_force (colvarvalue const &force)
 Apply the collective variable force, by communicating the atomic forces to the simulation program (Note: the member is not altered by this function).
virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const
 Square distance between x1 and x2 (can be redefined to transparently implement constraints, symmetries and periodicities).
virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const
 Gradient (with respect to x1) of the square distance (can be redefined to transparently implement constraints, symmetries and periodicities).
virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const
 Gradient (with respect to x2) of the square distance (can be redefined to transparently implement constraints, symmetries and periodicities).
virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const
 Return a positive number if x2>x1, zero if x2==x1, negative otherwise (can be redefined to transparently implement constraints, symmetries and periodicities) Note: it only works with scalar variables, otherwise raises an error.

Protected Attributes

cvm::atom_group atoms
 Atom group.
std::vector< cvm::atom_pos > ref_pos
 Reference coordinates.
cvm::rvector ref_pos_center
 Geometric center of the reference coordinates.
std::vector< cvm::rvector > eigenvec
 Eigenvector (of a normal or essential mode): will always have zero center.
cvm::real eigenvec_invnorm2
 Inverse square norm of the eigenvector.

Detailed Description

Colvar component: projection of 3N coordinates onto an eigenvector (colvarvalue::type_scalar type, range (-*:*)).

Definition at line 582 of file colvarcomp.h.

Constructor & Destructor Documentation

colvar::eigenvector::eigenvector (  std::string const &  conf  ) 
 

Constructor.

Definition at line 882 of file colvarcomp_distances.C.

References atoms, colvarmodule::atom_group::b_center, colvarmodule::atom_group::b_rotate, colvarmodule::atom_group::b_user_defined_fit, colvarmodule::rotation::calc_optimal_rotation(), colvarmodule::atom_group::center_ref_pos(), colvarmodule::atom_group::create_sorted_ids(), eigenvec, eigenvec_invnorm2, colvarmodule::fatal_error(), colvarmodule::load_coords(), colvarmodule::log(), colvarmodule::rvector::norm2(), colvar::cvc::parse_group(), colvarmodule::atom_group::ref_pos, ref_pos, ref_pos_center, colvarmodule::rotation::request_group1_gradients(), colvarmodule::rotation::request_group2_gradients(), colvarmodule::atom_group::rot, colvarmodule::rotation::rotate(), colvarmodule::atom_group::sorted_ids, and colvarvalue::type().

00883   : cvc (conf)
00884 {
00885   b_inverse_gradients = true;
00886   b_Jacobian_derivative = true;
00887   function_type = "eigenvector";
00888   x.type (colvarvalue::type_scalar);
00889 
00890   parse_group (conf, "atoms", atoms);
00891   atom_groups.push_back (&atoms);
00892 
00893 
00894   {
00895     bool const b_inline = get_keyval (conf, "refPositions", ref_pos, ref_pos);
00896 
00897     if (b_inline) {
00898       cvm::log ("Using reference positions from input file.\n");
00899       if (ref_pos.size() != atoms.size()) {
00900         cvm::fatal_error ("Error: reference positions do not "
00901                           "match the number of requested atoms.\n");
00902       }
00903     }
00904 
00905     std::string file_name;
00906     if (get_keyval (conf, "refPositionsFile", file_name)) {
00907 
00908       if (b_inline)
00909         cvm::fatal_error ("Error: refPositions and refPositionsFile cannot be specified at the same time.\n");
00910 
00911       std::string file_col;
00912       double file_col_value;
00913       if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) {
00914         // use PDB flags if column is provided
00915         bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0);
00916         if (found && !file_col_value)
00917           cvm::fatal_error ("Error: refPositionsColValue, "
00918                             "if provided, must be non-zero.\n");
00919       } else {
00920         // if not, use atom indices
00921         atoms.create_sorted_ids();
00922       }
00923 
00924       ref_pos.resize (atoms.size());
00925       cvm::load_coords (file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value);
00926     }
00927   }
00928 
00929   // save for later the geometric center of the provided positions (may not be the origin)
00930   cvm::rvector ref_pos_center (0.0, 0.0, 0.0);
00931   for (size_t i = 0; i < atoms.size(); i++) {
00932     ref_pos_center += ref_pos[i];
00933   }
00934   ref_pos_center *= 1.0 / atoms.size();
00935 
00936   if (atoms.b_user_defined_fit) {
00937     cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
00938   } else {
00939     // default: fit everything
00940     cvm::log ("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: "
00941               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
00942     atoms.b_center = true;
00943     atoms.b_rotate = true;
00944     atoms.ref_pos = ref_pos;
00945     atoms.center_ref_pos();
00946 
00947     // request the calculation of the derivatives of the rotation defined by the atom group
00948     atoms.rot.request_group1_gradients (atoms.size());
00949     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
00950     // this is only required for ABF, but we do both groups here for better caching
00951     atoms.rot.request_group2_gradients (atoms.size());
00952   }
00953 
00954   {
00955     bool const b_inline = get_keyval (conf, "vector", eigenvec, eigenvec);
00956     // now load the eigenvector
00957     if (b_inline) {
00958       cvm::log ("Using vector components from input file.\n");
00959       if (eigenvec.size() != atoms.size()) {
00960         cvm::fatal_error ("Error: vector components do not "
00961                           "match the number of requested atoms.\n");
00962       }
00963     }
00964 
00965     std::string file_name;
00966     if (get_keyval (conf, "vectorFile", file_name)) {
00967 
00968       if (b_inline)
00969         cvm::fatal_error ("Error: vector and vectorFile cannot be specified at the same time.\n");
00970 
00971       std::string file_col;
00972       double file_col_value;
00973       if (get_keyval (conf, "vectorCol", file_col, std::string (""))) {
00974         // use PDB flags if column is provided
00975         bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0);
00976         if (found && !file_col_value)
00977           cvm::fatal_error ("Error: vectorColValue, "
00978                             "if provided, must be non-zero.\n");
00979       } else {
00980         // if not, use atom indices
00981         atoms.create_sorted_ids();
00982       }
00983 
00984       eigenvec.resize (atoms.size());
00985       cvm::load_coords (file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value);
00986     }
00987   }
00988 
00989   if (!ref_pos.size() || !eigenvec.size()) {
00990     cvm::fatal_error ("Error: both reference coordinates"
00991                       "and eigenvector must be defined.\n");
00992   }
00993 
00994   cvm::rvector eig_center (0.0, 0.0, 0.0);
00995   for (size_t i = 0; i < atoms.size(); i++) {
00996     eig_center += eigenvec[i];
00997   }
00998   eig_center *= 1.0 / atoms.size();
00999   cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n");
01000 
01001   bool b_difference_vector = false;
01002   get_keyval (conf, "differenceVector", b_difference_vector, false);
01003 
01004   if (b_difference_vector) {
01005 
01006     if (atoms.b_center) {
01007       // both sets should be centered on the origin for fitting
01008       for (size_t i = 0; i < atoms.size(); i++) {
01009         eigenvec[i] -= eig_center;
01010         ref_pos[i]  -= ref_pos_center;
01011       }
01012     }
01013     if (atoms.b_rotate) {
01014       atoms.rot.calc_optimal_rotation (eigenvec, ref_pos);
01015       for (size_t i = 0; i < atoms.size(); i++) {
01016         eigenvec[i] = atoms.rot.rotate (eigenvec[i]);
01017       }
01018     }
01019     cvm::log ("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n");
01020     for (size_t i = 0; i < atoms.size(); i++) {
01021       eigenvec[i] -= ref_pos[i];
01022     }
01023     if (atoms.b_center) {
01024       // bring back the ref positions to where they were
01025       for (size_t i = 0; i < atoms.size(); i++) {
01026         ref_pos[i] += ref_pos_center;
01027       }
01028     }
01029 
01030   } else {
01031     cvm::log ("Centering the provided vector to zero.\n");
01032     for (size_t i = 0; i < atoms.size(); i++) {
01033       eigenvec[i] -= eig_center;
01034     }
01035   }
01036 
01037   // cvm::log ("The first three components (v1x, v1y, v1z) of the resulting vector are: "+cvm::to_str (eigenvec[0])+".\n");
01038 
01039   // for inverse gradients
01040   eigenvec_invnorm2 = 0.0;
01041   for (size_t i = 0; i < atoms.size(); i++) {
01042     eigenvec_invnorm2 += eigenvec[i].norm2();
01043   }
01044   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
01045 
01046   if (b_difference_vector) {
01047     cvm::log ("\"differenceVector\" is on: normalizing the vector.\n");
01048     for (size_t i = 0; i < atoms.size(); i++) {
01049       eigenvec[i] *= eigenvec_invnorm2;
01050     }
01051   } else {
01052     cvm::log ("The norm of the vector is |v| = "+cvm::to_str (eigenvec_invnorm2)+".\n");
01053   }
01054 }

virtual colvar::eigenvector::~eigenvector (   )  [inline, virtual]
 

Definition at line 606 of file colvarcomp.h.

00606 {}

Member Function Documentation

void colvar::eigenvector::apply_force (  colvarvalue const &  force  )  [virtual]
 

Apply the collective variable force, by communicating the atomic forces to the simulation program (Note: the member is not altered by this function).

Note: multiple calls to this function within the same simulation step will add the forces altogether

Parameters:
cvforce The collective variable force, usually coming from the biases and eventually manipulated by the parent object

Implements colvar::cvc.

Definition at line 1079 of file colvarcomp_distances.C.

References colvarmodule::atom_group::apply_colvar_force(), atoms, colvarmodule::atom_group::noforce, and colvarvalue::real_value.

01080 {
01081   if (!atoms.noforce)
01082     atoms.apply_colvar_force (force.real_value);
01083 }

void colvar::eigenvector::calc_force_invgrads (   )  [virtual]
 

Calculate the total force from the system using the inverse atomic gradients.

Reimplemented from colvar::cvc.

Definition at line 1086 of file colvarcomp_distances.C.

References atoms, eigenvec_invnorm2, colvarmodule::atom_group::read_system_forces(), colvarvalue::real_value, and colvarmodule::atom_group::system_force().

01087 {
01088   atoms.read_system_forces();
01089   ft.real_value = 0.0;
01090     
01091   for (size_t ia = 0; ia < atoms.size(); ia++) {
01092     ft.real_value += eigenvec_invnorm2 * atoms[ia].grad *
01093       atoms[ia].system_force;
01094   }
01095 }

void colvar::eigenvector::calc_gradients (   )  [virtual]
 

Calculate the atomic gradients, to be reused later in order to apply forces.

Implements colvar::cvc.

Definition at line 1066 of file colvarcomp_distances.C.

References atoms, colvar::cvc::debug_gradients(), eigenvec, and colvarmodule::log().

01067 {
01068   for (size_t ia = 0; ia < atoms.size(); ia++) {
01069     atoms[ia].grad = eigenvec[ia];
01070   }
01071 
01072   if (b_debug_gradients) {
01073     cvm::log ("Debugging gradients:\n");
01074     debug_gradients (atoms);
01075   }
01076 }

void colvar::eigenvector::calc_Jacobian_derivative (   )  [virtual]
 

Calculate the divergence of the inverse atomic gradients.

Reimplemented from colvar::cvc.

Definition at line 1098 of file colvarcomp_distances.C.

References colvarmodule::atom_pos, atoms, colvarmodule::rotation::dQ0_1, eigenvec, eigenvec_invnorm2, j, colvarmodule::rotation::q, colvarvalue::real_value, and colvarmodule::atom_group::rot.

01099 {
01100   // gradient of the rotation matrix
01101   cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
01102   cvm::quaternion &quat0 = atoms.rot.q;
01103 
01104   // gradients of products of 2 quaternion components 
01105   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
01106 
01107   cvm::atom_pos x_relative; 
01108   cvm::real sum = 0.0;
01109 
01110   for (size_t ia = 0; ia < atoms.size(); ia++) {
01111 
01112     // Gradient of optimal quaternion wrt current Cartesian position
01113     // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t
01114     // we can just transpose the derivatives of the direct matrix
01115     cvm::vector1d< cvm::rvector, 4 >      &dq_1 = atoms.rot.dQ0_1[ia];
01116 
01117     g11 = 2.0 * quat0[1]*dq_1[1];
01118     g22 = 2.0 * quat0[2]*dq_1[2];
01119     g33 = 2.0 * quat0[3]*dq_1[3];
01120     g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0];
01121     g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0];
01122     g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0];
01123     g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1];
01124     g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1];
01125     g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2];
01126 
01127     // Gradient of the inverse rotation matrix wrt current Cartesian position
01128     // (transpose of the gradient of the direct rotation)
01129     grad_rot_mat[0][0] = -2.0 * (g22 + g33); 
01130     grad_rot_mat[0][1] =  2.0 * (g12 + g03); 
01131     grad_rot_mat[0][2] =  2.0 * (g13 - g02); 
01132     grad_rot_mat[1][0] =  2.0 * (g12 - g03); 
01133     grad_rot_mat[1][1] = -2.0 * (g11 + g33); 
01134     grad_rot_mat[1][2] =  2.0 * (g01 + g23); 
01135     grad_rot_mat[2][0] =  2.0 * (g02 + g13); 
01136     grad_rot_mat[2][1] =  2.0 * (g23 - g01); 
01137     grad_rot_mat[2][2] = -2.0 * (g11 + g22); 
01138 
01139     for (size_t i = 0; i < 3; i++) {
01140       for (size_t j = 0; j < 3; j++) {
01141         sum += grad_rot_mat[i][j][i] * eigenvec[ia][j];
01142       }
01143     }
01144   }
01145 
01146   jd.real_value = sum * std::sqrt (eigenvec_invnorm2); 
01147 }

void colvar::eigenvector::calc_value (   )  [virtual]
 

Calculate the variable.

Implements colvar::cvc.

Definition at line 1057 of file colvarcomp_distances.C.

References atoms, eigenvec, colvarvalue::real_value, and ref_pos.

01058 {
01059   x.real_value = 0.0;
01060   for (size_t i = 0; i < atoms.size(); i++) {
01061     x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i];
01062   }
01063 }

virtual cvm::real colvar::eigenvector::compare (  colvarvalue const &  x1,
colvarvalue const &  x2
)  const [virtual]
 

Return a positive number if x2>x1, zero if x2==x1, negative otherwise (can be redefined to transparently implement constraints, symmetries and periodicities) Note: it only works with scalar variables, otherwise raises an error.

Reimplemented from colvar::cvc.

virtual cvm::real colvar::eigenvector::dist2 (  colvarvalue const &  x1,
colvarvalue const &  x2
)  const [virtual]
 

Square distance between x1 and x2 (can be redefined to transparently implement constraints, symmetries and periodicities).

colvar::cvc::dist2() and the related functions are declared as "const" functions, but not "static", because additional parameters defining the metrics (e.g. the periodicity) may be specific to each colvar::cvc object.

If symmetries or periodicities are present, the colvar::cvc::dist2() should be redefined to return the "closest distance" value and colvar::cvc::dist2_lgrad(), colvar::cvc::dist2_rgrad() to return its gradients.

If constraints are present (and not already implemented by any of the types), the colvar::cvc::dist2_lgrad() and colvar::cvc::dist2_rgrad() functions should be redefined to provide a gradient which is compatible with the constraint, i.e. already deprived of its component normal to the constraint hypersurface.

Finally, another useful application, if you are performing very many operations with these functions, could be to override the member functions and access directly its member data. For instance: to define dist2(x1,x2) as (x2.real_value-x1.real_value)*(x2.real_value-x1.real_value) in case of a scalar type.

Reimplemented from colvar::cvc.

virtual colvarvalue colvar::eigenvector::dist2_lgrad (  colvarvalue const &  x1,
colvarvalue const &  x2
)  const [virtual]
 

Gradient (with respect to x1) of the square distance (can be redefined to transparently implement constraints, symmetries and periodicities).

Reimplemented from colvar::cvc.

virtual colvarvalue colvar::eigenvector::dist2_rgrad (  colvarvalue const &  x1,
colvarvalue const &  x2
)  const [virtual]
 

Gradient (with respect to x2) of the square distance (can be redefined to transparently implement constraints, symmetries and periodicities).

Reimplemented from colvar::cvc.

Member Data Documentation

cvm::atom_group colvar::eigenvector::atoms [protected]
 

Atom group.

Definition at line 588 of file colvarcomp.h.

Referenced by apply_force(), calc_force_invgrads(), calc_gradients(), calc_Jacobian_derivative(), calc_value(), and eigenvector().

std::vector<cvm::rvector> colvar::eigenvector::eigenvec [protected]
 

Eigenvector (of a normal or essential mode): will always have zero center.

Definition at line 597 of file colvarcomp.h.

Referenced by calc_gradients(), calc_Jacobian_derivative(), calc_value(), and eigenvector().

cvm::real colvar::eigenvector::eigenvec_invnorm2 [protected]
 

Inverse square norm of the eigenvector.

Definition at line 600 of file colvarcomp.h.

Referenced by calc_force_invgrads(), calc_Jacobian_derivative(), and eigenvector().

std::vector<cvm::atom_pos> colvar::eigenvector::ref_pos [protected]
 

Reference coordinates.

Definition at line 591 of file colvarcomp.h.

Referenced by calc_value(), and eigenvector().

cvm::rvector colvar::eigenvector::ref_pos_center [protected]
 

Geometric center of the reference coordinates.

Definition at line 594 of file colvarcomp.h.

Referenced by eigenvector().

The documentation for this class was generated from the following files:
Generated on Tue May 21 04:07:24 2013 for NAMD by  doxygen 1.3.9.1