Difference for src/colvarcomp_rotations.C from version 1.15 to 1.16

version 1.15version 1.16
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include <cmath> #include <cmath>
  
Line 14
Line 14
   : 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");
     }     }
Line 42
Line 41
                             "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);
     }     }
   }   }
  
Line 58
Line 57
   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());
   }   }
  
 } }
Line 88
Line 87
  
 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;
Line 116
Line 113
 { {
   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]);
       }       }
     }     }
   }   }
Line 145
Line 142
  
 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);
Line 167
Line 161
       ((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); 
   }   }
 } }
  
Line 180
Line 170
 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);
   }   }
 } }
  
Line 205
Line 195
  
 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;
 } }
  
Line 214
Line 204
 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); 
   }   }
 } }
  
Line 228
Line 214
 { {
   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);
   }   }
 } }
  
Line 261
Line 247
  
 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);
 } }
Line 276
Line 259
 { {
   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); 
   }   }
 } }
  
Line 294
Line 272
 { {
   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);
   }   }
 } }
  
Line 331
Line 309
  
 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);
Line 347
Line 322
 { {
   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]);
     }     }
   }   }
 } }
Line 360
Line 335
 { {
   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);
   }   }
 } }

Legend:
Removed in v.1.15 
changed lines
 Added in v.1.16


Made by using version 1.53 of cvs2html