Difference for src/colvarcomp_angles.C from version 1.12 to 1.13

version 1.12version 1.13
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include "colvarmodule.h" #include "colvarmodule.h"
 #include "colvar.h" #include "colvar.h"
Line 11
Line 11
   : cvc(conf)   : cvc(conf)
 { {
   function_type = "angle";   function_type = "angle";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
   parse_group(conf, "group1", group1);   provide(f_cvc_com_based);
   parse_group(conf, "group2", group2); 
   parse_group(conf, "group3", group3);   group1 = parse_group(conf, "group1");
   atom_groups.push_back(&group1);   group2 = parse_group(conf, "group2");
   atom_groups.push_back(&group2);   group3 = parse_group(conf, "group3");
   atom_groups.push_back(&group3); 
   if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {   init_total_force_params(conf);
     cvm::log("Computing system force on group 1 only"); 
   } 
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 29
Line 28
 colvar::angle::angle(cvm::atom const &a1, colvar::angle::angle(cvm::atom const &a1,
                       cvm::atom const &a2,                       cvm::atom const &a2,
                       cvm::atom const &a3)                       cvm::atom const &a3)
   : group1(std::vector<cvm::atom> (1, a1)), 
     group2(std::vector<cvm::atom> (1, a2)), 
     group3(std::vector<cvm::atom> (1, a3)) 
 { {
   function_type = "angle";   function_type = "angle";
   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);
   atom_groups.push_back(&group1); 
   atom_groups.push_back(&group2);   group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
   atom_groups.push_back(&group3);   group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
    group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
    atom_groups.push_back(group1);
    atom_groups.push_back(group2);
    atom_groups.push_back(group3);
  
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
Line 54
Line 54
  
 void colvar::angle::calc_value() void colvar::angle::calc_value()
 { {
   cvm::atom_pos const g1_pos = group1.center_of_mass();   cvm::atom_pos const g1_pos = group1->center_of_mass();
   cvm::atom_pos const g2_pos = group2.center_of_mass();   cvm::atom_pos const g2_pos = group2->center_of_mass();
   cvm::atom_pos const g3_pos = group3.center_of_mass();   cvm::atom_pos const g3_pos = group3->center_of_mass();
  
   r21  = cvm::position_distance(g2_pos, g1_pos);   r21  = cvm::position_distance(g2_pos, g1_pos);
   r21l = r21.norm();   r21l = r21.norm();
Line 71
Line 71
  
 void colvar::angle::calc_gradients() void colvar::angle::calc_gradients()
 { {
   size_t i; 
   cvm::real const cos_theta = (r21*r23)/(r21l*r23l);   cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
   cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta);   cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta);
  
Line 81
Line 80
   dxdr3 = (180.0/PI) * dxdcos *   dxdr3 = (180.0/PI) * dxdcos *
     (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );     (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
  
   for (i = 0; i < group1.size(); i++) {   group1->set_weighted_gradient(dxdr1);
     group1[i].grad = (group1[i].mass/group1.total_mass) *   group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0));
       (dxdr1);   group3->set_weighted_gradient(dxdr3);
   } 
  
   for (i = 0; i < group2.size(); i++) { 
     group2[i].grad = (group2[i].mass/group2.total_mass) * 
       (dxdr1 + dxdr3) * (-1.0); 
   } 
  
   for (i = 0; i < group3.size(); i++) { 
     group3[i].grad = (group3[i].mass/group3.total_mass) * 
       (dxdr3); 
   } 
 } }
  
 void colvar::angle::calc_force_invgrads() void colvar::angle::calc_force_invgrads()
Line 105
Line 93
   // centered on group2, which means group2 is kept fixed   // centered on group2, which means group2 is kept fixed
   // when propagating changes in the angle)   // when propagating changes in the angle)
  
   if (b_1site_force) {   if (is_enabled(f_cvc_one_site_total_force)) {
     group1.read_system_forces();     group1->read_total_forces();
     cvm::real norm_fact = 1.0 / dxdr1.norm2();     cvm::real norm_fact = 1.0 / dxdr1.norm2();
     ft.real_value = norm_fact * dxdr1 * group1.system_force();     ft.real_value = norm_fact * dxdr1 * group1->total_force();
   } else {   } else {
     group1.read_system_forces();     group1->read_total_forces();
     group3.read_system_forces();     group3->read_total_forces();
     cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2());     cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2());
     ft.real_value = norm_fact * ( dxdr1 * group1.system_force()     ft.real_value = norm_fact * ( dxdr1 * group1->total_force()
                                 + dxdr3 * group3.system_force());                                 + dxdr3 * group3->total_force());
   }   }
   return;   return;
 } }
Line 130
Line 118
  
 void colvar::angle::apply_force(colvarvalue const &force) void colvar::angle::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)
      group2->apply_colvar_force(force.real_value);
  
    if (!group3->noforce)
      group3->apply_colvar_force(force.real_value);
  }
  
  
  
  
  colvar::dipole_angle::dipole_angle(std::string const &conf)
    : cvc(conf)
  {
    function_type = "dipole_angle";
    group1 = parse_group(conf, "group1");
    group2 = parse_group(conf, "group2");
    group3 = parse_group(conf, "group3");
  
    init_total_force_params(conf);
  
    x.type(colvarvalue::type_scalar);
  }
  
  
  colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
                        cvm::atom const &a2,
                        cvm::atom const &a3)
  {
    function_type = "dipole_angle";
  
    group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
    group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
    group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
    atom_groups.push_back(group1);
    atom_groups.push_back(group2);
    atom_groups.push_back(group3);
  
    x.type(colvarvalue::type_scalar);
  }
  
   if (!group2.noforce) 
     group2.apply_colvar_force(force.real_value); 
  
   if (!group3.noforce) colvar::dipole_angle::dipole_angle()
     group3.apply_colvar_force(force.real_value); {
    function_type = "dipole_angle";
    x.type(colvarvalue::type_scalar);
  }
  
  
  void colvar::dipole_angle::calc_value()
  {
    cvm::atom_pos const g1_pos = group1->center_of_mass();
    cvm::atom_pos const g2_pos = group2->center_of_mass();
    cvm::atom_pos const g3_pos = group3->center_of_mass();
  
    group1->calc_dipole(g1_pos);
  
    r21 = group1->dipole();
    r21l = r21.norm();
    r23  = cvm::position_distance(g2_pos, g3_pos);
    r23l = r23.norm();
  
    cvm::real     const cos_theta = (r21*r23)/(r21l*r23l);
  
    x.real_value = (180.0/PI) * std::acos(cos_theta);
  }
  
  //to be implemented
  //void colvar::dipole_angle::calc_force_invgrads(){}
  //void colvar::dipole_angle::calc_Jacobian_derivative(){}
  
  void colvar::dipole_angle::calc_gradients()
  {
    cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
    cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta);
  
    dxdr1 = (180.0/PI) * dxdcos *
    (1.0/r21l)* (r23/r23l + (-1.0) * cos_theta * r21/r21l );
  
    dxdr3 =  (180.0/PI) * dxdcos *
      (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
  
    //this auxiliar variables are to avoid numerical errors inside "for"
    double aux1 = group1->total_charge/group1->total_mass;
    // double aux2 = group2->total_charge/group2->total_mass;
    // double aux3 = group3->total_charge/group3->total_mass;
  
    size_t i;
    for (i = 0; i < group1->size(); i++) {
      (*group1)[i].grad =((*group1)[i].charge + (-1)* (*group1)[i].mass * aux1) * (dxdr1);
    }
  
    for (i = 0; i < group2->size(); i++) {
      (*group2)[i].grad = ((*group2)[i].mass/group2->total_mass)* dxdr3 * (-1.0);
    }
  
    for (i = 0; i < group3->size(); i++) {
      (*group3)[i].grad =((*group3)[i].mass/group3->total_mass) * (dxdr3);
    }
  }
  
  
  void colvar::dipole_angle::apply_force(colvarvalue const &force)
  {
    if (!group1->noforce)
      group1->apply_colvar_force(force.real_value);
  
    if (!group2->noforce)
      group2->apply_colvar_force(force.real_value);
  
    if (!group3->noforce)
      group3->apply_colvar_force(force.real_value);
 } }
  
  
Line 149
Line 243
   function_type = "dihedral";   function_type = "dihedral";
   period = 360.0;   period = 360.0;
   b_periodic = true;   b_periodic = true;
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
   if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {   provide(f_cvc_com_based);
     cvm::log("Computing system force on group 1 only"); 
   }   group1 = parse_group(conf, "group1");
   parse_group(conf, "group1", group1);   group2 = parse_group(conf, "group2");
   parse_group(conf, "group2", group2);   group3 = parse_group(conf, "group3");
   parse_group(conf, "group3", group3);   group4 = parse_group(conf, "group4");
   parse_group(conf, "group4", group4); 
   atom_groups.push_back(&group1);   init_total_force_params(conf);
   atom_groups.push_back(&group2); 
   atom_groups.push_back(&group3); 
   atom_groups.push_back(&group4); 
  
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
Line 171
Line 262
                             cvm::atom const &a2,                             cvm::atom const &a2,
                             cvm::atom const &a3,                             cvm::atom const &a3,
                             cvm::atom const &a4)                             cvm::atom const &a4)
   : group1(std::vector<cvm::atom> (1, a1)), 
     group2(std::vector<cvm::atom> (1, a2)), 
     group3(std::vector<cvm::atom> (1, a3)), 
     group4(std::vector<cvm::atom> (1, a4)) 
 { {
   if (cvm::debug())   if (cvm::debug())
     cvm::log("Initializing dihedral object from atom groups.\n");     cvm::log("Initializing dihedral object from atom groups.\n");
Line 182
Line 269
   function_type = "dihedral";   function_type = "dihedral";
   period = 360.0;   period = 360.0;
   b_periodic = true;   b_periodic = true;
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
    provide(f_cvc_com_based);
  
   b_1site_force = false;   b_1site_force = false;
  
   atom_groups.push_back(&group1);   group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
   atom_groups.push_back(&group2);   group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
   atom_groups.push_back(&group3);   group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
   atom_groups.push_back(&group4);   group4 = new cvm::atom_group(std::vector<cvm::atom>(1, a4));
    atom_groups.push_back(group1);
    atom_groups.push_back(group2);
    atom_groups.push_back(group3);
    atom_groups.push_back(group4);
  
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
  
Line 203
Line 296
   function_type = "dihedral";   function_type = "dihedral";
   period = 360.0;   period = 360.0;
   b_periodic = true;   b_periodic = true;
   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::dihedral::calc_value() void colvar::dihedral::calc_value()
 { {
   cvm::atom_pos const g1_pos = group1.center_of_mass();   cvm::atom_pos const g1_pos = group1->center_of_mass();
   cvm::atom_pos const g2_pos = group2.center_of_mass();   cvm::atom_pos const g2_pos = group2->center_of_mass();
   cvm::atom_pos const g3_pos = group3.center_of_mass();   cvm::atom_pos const g3_pos = group3->center_of_mass();
   cvm::atom_pos const g4_pos = group4.center_of_mass();   cvm::atom_pos const g4_pos = group4->center_of_mass();
  
   // Usual sign convention: r12 = r2 - r1   // Usual sign convention: r12 = r2 - r1
   r12 = cvm::position_distance(g1_pos, g2_pos);   r12 = cvm::position_distance(g1_pos, g2_pos);
Line 299
Line 392
               +dsindB.y*r34.x - dsindB.x*r34.y);               +dsindB.y*r34.x - dsindB.x*r34.y);
   }   }
  
   size_t i;   group1->set_weighted_gradient(-f1);
   for (i = 0; i < group1.size(); i++)   group2->set_weighted_gradient(-f2 + f1);
     group1[i].grad = (group1[i].mass/group1.total_mass) * (-f1);   group3->set_weighted_gradient(-f3 + f2);
    group4->set_weighted_gradient(f3);
   for (i = 0; i < group2.size(); i++) 
     group2[i].grad = (group2[i].mass/group2.total_mass) * (-f2 + f1); 
  
   for (i = 0; i < group3.size(); i++) 
  group3[i].grad = (group3[i].mass/group3.total_mass) * (-f3 + f2); 
  
   for (i = 0; i < group4.size(); i++) 
     group4[i].grad = (group4[i].mass/group4.total_mass) * (f3); 
 } }
  
  
Line 332
Line 417
   cvm::real const fact1 = d12 * std::sqrt(1.0 - dot1 * dot1);   cvm::real const fact1 = d12 * std::sqrt(1.0 - dot1 * dot1);
   cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4);   cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4);
  
   group1.read_system_forces();   group1->read_total_forces();
   if ( b_1site_force ) {   if (is_enabled(f_cvc_one_site_total_force)) {
     // This is only measuring the force on group 1     // This is only measuring the force on group 1
     ft.real_value = PI/180.0 * fact1 * (cross1 * group1.system_force());     ft.real_value = PI/180.0 * fact1 * (cross1 * group1->total_force());
   } else {   } else {
     // Default case: use groups 1 and 4     // Default case: use groups 1 and 4
     group4.read_system_forces();     group4->read_total_forces();
     ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1.system_force())     ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1->total_force())
        + fact4 * (cross4 * group4.system_force()));        + fact4 * (cross4 * group4->total_force()));
   }   }
 } }
  
Line 354
Line 439
  
 void colvar::dihedral::apply_force(colvarvalue const &force) void colvar::dihedral::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);
  
   if (!group3.noforce)   if (!group3->noforce)
     group3.apply_colvar_force(force.real_value);     group3->apply_colvar_force(force.real_value);
  
   if (!group4.noforce)   if (!group4->noforce)
     group4.apply_colvar_force(force.real_value);     group4->apply_colvar_force(force.real_value);
 } }
  
  

Legend:
Removed in v.1.12 
changed lines
 Added in v.1.13


Made by using version 1.53 of cvs2html