Difference for src/colvarcomp_coordnums.C from version 1.11 to 1.12

version 1.11version 1.12
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include <cmath> #include <cmath>
  
Line 71
Line 71
 } }
  
  
  
 colvar::coordnum::coordnum(std::string const &conf) colvar::coordnum::coordnum(std::string const &conf)
   : distance(conf), b_anisotropic(false), b_group2_center_only(false)   : cvc(conf), b_anisotropic(false), b_group2_center_only(false)
 { {
   function_type = "coordnum";   function_type = "coordnum";
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
  
   // group1 and group2 are already initialized by distance()   group1 = parse_group(conf, "group1");
   if (group1.b_dummy)   group2 = parse_group(conf, "group2");
     cvm::fatal_error("Error: only group2 is allowed to be a dummy atom\n"); 
  
    if (group1->b_dummy)
      cvm::fatal_error("Error: only group2 is allowed to be a dummy atom\n");
  
   // need to specify this explicitly because the distance() constructor   bool const b_isotropic = get_keyval(conf, "cutoff", r0,
   // has set it to true 
   b_inverse_gradients = false; 
  
   bool const b_scale = get_keyval(conf, "cutoff", r0, 
                                    cvm::real(4.0 * cvm::unit_angstrom()));                                    cvm::real(4.0 * cvm::unit_angstrom()));
  
   if (get_keyval(conf, "cutoff3", r0_vec,   if (get_keyval(conf, "cutoff3", r0_vec, cvm::rvector(4.0 * cvm::unit_angstrom(),
                   cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {                                                        4.0 * cvm::unit_angstrom(),
                                                         4.0 * cvm::unit_angstrom()))) {
      if (b_isotropic) {
        cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
                   INPUT_ERROR);
      }
  
     if (b_scale) 
       cvm::fatal_error("Error: cannot specify \"scale\" and " 
                         "\"scale3\" at the same time.\n"); 
     b_anisotropic = true;     b_anisotropic = true;
     // remove meaningless negative signs     // remove meaningless negative signs
     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
Line 107
Line 105
   get_keyval(conf, "expDenom", ed, int(12));   get_keyval(conf, "expDenom", ed, int(12));
  
   if ( (en%2) || (ed%2) ) {   if ( (en%2) || (ed%2) ) {
     cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");     cvm::error("Error: odd exponents provided, can only use even ones.\n", INPUT_ERROR);
   }   }
  
   get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2.b_dummy);   get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
 } }
  
  
Line 130
Line 128
  
     // create a fake atom to hold the group2 com coordinates     // create a fake atom to hold the group2 com coordinates
     cvm::atom group2_com_atom;     cvm::atom group2_com_atom;
     group2_com_atom.pos = group2.center_of_mass();     group2_com_atom.pos = group2->center_of_mass();
  
     if (b_anisotropic) {     if (b_anisotropic) {
       for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
         x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);         x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);
     } else {     } else {
       for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
         x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);         x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);
     }     }
  
   } else {   } else {
  
     if (b_anisotropic) {     if (b_anisotropic) {
       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++) {
           x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);           x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);
         }         }
     } 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++) {
           x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);           x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);
         }         }
     }     }
Line 163
Line 161
  
     // create a fake atom to hold the group2 com coordinates     // create a fake atom to hold the group2 com coordinates
     cvm::atom group2_com_atom;     cvm::atom group2_com_atom;
     group2_com_atom.pos = group2.center_of_mass();     group2_com_atom.pos = group2->center_of_mass();
  
  
     if (b_anisotropic) {     if (b_anisotropic) {
       for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
         switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);         switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);
     } else {     } else {
       for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
         switching_function<true>(r0, en, ed, *ai1, group2_com_atom);         switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
     }     }
  
     group2.set_weighted_gradient(group2_com_atom.grad);     group2->set_weighted_gradient(group2_com_atom.grad);
  
   } else {   } else {
  
     if (b_anisotropic) {     if (b_anisotropic) {
       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++) {
           switching_function<true>(r0_vec, en, ed, *ai1, *ai2);           switching_function<true>(r0_vec, en, ed, *ai1, *ai2);
         }         }
     } 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++) {
           switching_function<true>(r0, en, ed, *ai1, *ai2);           switching_function<true>(r0, en, ed, *ai1, *ai2);
         }         }
     }     }
Line 194
Line 192
  
 void colvar::coordnum::apply_force(colvarvalue const &force) void colvar::coordnum::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);
 } }
  
  
Line 289
Line 287
  
  
 colvar::selfcoordnum::selfcoordnum(std::string const &conf) colvar::selfcoordnum::selfcoordnum(std::string const &conf)
  : distance(conf, false)   : cvc(conf)
 { {
   function_type = "selfcoordnum";   function_type = "selfcoordnum";
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
  
   // group1 is already initialized by distance()   group1 = parse_group(conf, "group1");
  
   // need to specify this explicitly because the distance() constructor 
   // has set it to true 
   b_inverse_gradients = false; 
  
   get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom()));   get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom()));
   get_keyval(conf, "expNumer", en, int(6) );   get_keyval(conf, "expNumer", en, int(6) );
Line 320
Line 314
 void colvar::selfcoordnum::calc_value() void colvar::selfcoordnum::calc_value()
 { {
   x.real_value = 0.0;   x.real_value = 0.0;
   for (size_t i = 0; i < group1.size() - 1; i++) {   for (size_t i = 0; i < group1->size() - 1; i++) {
     for (size_t j = i + 1; j < group1.size(); j++) {     for (size_t j = i + 1; j < group1->size(); j++) {
       x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, group1[i], group1[j]);       x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, (*group1)[i], (*group1)[j]);
     }     }
   }   }
 } }
Line 330
Line 324
  
 void colvar::selfcoordnum::calc_gradients() void colvar::selfcoordnum::calc_gradients()
 { {
   for (size_t i = 0; i < group1.size() - 1; i++) {   for (size_t i = 0; i < group1->size() - 1; i++) {
     for (size_t j = i + 1; j < group1.size(); j++) {     for (size_t j = i + 1; j < group1->size(); j++) {
       colvar::coordnum::switching_function<true>(r0, en, ed, group1[i], group1[j]);       colvar::coordnum::switching_function<true>(r0, en, ed, (*group1)[i], (*group1)[j]);
     }     }
   }   }
 } }
  
 void colvar::selfcoordnum::apply_force(colvarvalue const &force) void colvar::selfcoordnum::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);
   }   }
 } }
  
  
  // groupcoordnum member functions
  colvar::groupcoordnum::groupcoordnum(std::string const &conf)
    : distance(conf), b_anisotropic(false)
  {
    function_type = "groupcoordnum";
    x.type(colvarvalue::type_scalar);
  
    // group1 and group2 are already initialized by distance()
    if (group1->b_dummy || group2->b_dummy)
      cvm::fatal_error("Error: neither group can be a dummy atom\n");
  
    bool const b_scale = get_keyval(conf, "cutoff", r0,
                                    cvm::real(4.0 * cvm::unit_angstrom()));
  
    if (get_keyval(conf, "cutoff3", r0_vec,
                   cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
  
      if (b_scale)
        cvm::fatal_error("Error: cannot specify \"scale\" and "
                         "\"scale3\" at the same time.\n");
      b_anisotropic = true;
      // remove meaningless negative signs
      if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
      if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
      if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
    }
  
    get_keyval(conf, "expNumer", en, int(6) );
    get_keyval(conf, "expDenom", ed, int(12));
  
    if ( (en%2) || (ed%2) ) {
      cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
    }
  
  }
  
  
  colvar::groupcoordnum::groupcoordnum()
    : b_anisotropic(false)
  {
    function_type = "groupcoordnum";
    x.type(colvarvalue::type_scalar);
  }
  
  
  template<bool calculate_gradients>
  cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
                                                      int const &en,
                                                      int const &ed,
                                                      cvm::atom &A1,
                                                      cvm::atom &A2)
  {
    cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
    cvm::real const l2 = diff.norm2()/(r0*r0);
  
    // Assume en and ed are even integers, and avoid sqrt in the following
    int const en2 = en/2;
    int const ed2 = ed/2;
  
    cvm::real const xn = std::pow(l2, en2);
    cvm::real const xd = std::pow(l2, ed2);
    cvm::real const func = (1.0-xn)/(1.0-xd);
  
    if (calculate_gradients) {
      cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
      cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
      A1.grad += (-1.0)*dFdl2*dl2dx;
      A2.grad +=        dFdl2*dl2dx;
    }
  
    return func;
  }
  
  
  #if 0  // AMG: I don't think there's any reason to support anisotropic,
         //      and I don't have those flags below in calc_value, but
         //      if I need them, I'll also need to uncomment this method
  template<bool calculate_gradients>
  cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
                                                      int const &en,
                                                      int const &ed,
                                                      cvm::atom &A1,
                                                      cvm::atom &A2)
  {
    cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
    cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
    cvm::real const l2 = scal_diff.norm2();
  
    // Assume en and ed are even integers, and avoid sqrt in the following
    int const en2 = en/2;
    int const ed2 = ed/2;
  
    cvm::real const xn = std::pow(l2, en2);
    cvm::real const xd = std::pow(l2, ed2);
    cvm::real const func = (1.0-xn)/(1.0-xd);
  
    if (calculate_gradients) {
      cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
      cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
                               (2.0/(r0_vec.y*r0_vec.y))*diff.y,
                               (2.0/(r0_vec.z*r0_vec.z))*diff.z);
      A1.grad += (-1.0)*dFdl2*dl2dx;
      A2.grad +=        dFdl2*dl2dx;
    }
    return func;
  }
  #endif
  
  
  
  void colvar::groupcoordnum::calc_value()
  {
  
    // create fake atoms to hold the com coordinates
    cvm::atom group1_com_atom;
    cvm::atom group2_com_atom;
    group1_com_atom.pos = group1->center_of_mass();
    group2_com_atom.pos = group2->center_of_mass();
  
    x.real_value = coordnum::switching_function<false>(r0, en, ed,
                                                       group1_com_atom, group2_com_atom);
  
  }
  
  
  void colvar::groupcoordnum::calc_gradients()
  {
    cvm::atom group1_com_atom;
    cvm::atom group2_com_atom;
    group1_com_atom.pos = group1->center_of_mass();
    group2_com_atom.pos = group2->center_of_mass();
  
    coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
    group1->set_weighted_gradient(group1_com_atom.grad);
    group2->set_weighted_gradient(group2_com_atom.grad);
  
  }
  
  
  void colvar::groupcoordnum::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);
  }

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


Made by using version 1.53 of cvs2html