00001 #if (__cplusplus >= 201103L) 00002 00003 // This file is part of the Collective Variables module (Colvars). 00004 // The original version of Colvars and its updates are located at: 00005 // https://github.com/Colvars/colvars 00006 // Please update all Colvars source files before making any changes. 00007 // If you wish to distribute your changes, please submit them to the 00008 // Colvars repository at GitHub. 00009 00010 #include "colvarcomp.h" 00011 00012 colvar::linearCombination::linearCombination(std::string const &conf): cvc(conf) { 00013 // Lookup all available sub-cvcs 00014 for (auto it_cv_map = colvar::get_global_cvc_map().begin(); it_cv_map != colvar::get_global_cvc_map().end(); ++it_cv_map) { 00015 if (key_lookup(conf, it_cv_map->first.c_str())) { 00016 std::vector<std::string> sub_cvc_confs; 00017 get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs); 00018 for (auto it_sub_cvc_conf = sub_cvc_confs.begin(); it_sub_cvc_conf != sub_cvc_confs.end(); ++it_sub_cvc_conf) { 00019 cv.push_back((it_cv_map->second)(*(it_sub_cvc_conf))); 00020 } 00021 } 00022 } 00023 // Sort all sub CVs by their names 00024 std::sort(cv.begin(), cv.end(), colvar::compare_cvc); 00025 for (auto it_sub_cv = cv.begin(); it_sub_cv != cv.end(); ++it_sub_cv) { 00026 for (auto it_atom_group = (*it_sub_cv)->atom_groups.begin(); it_atom_group != (*it_sub_cv)->atom_groups.end(); ++it_atom_group) { 00027 register_atom_group(*it_atom_group); 00028 } 00029 } 00030 // Show useful error messages and prevent crashes if no sub CVC is found 00031 if (cv.size() == 0) { 00032 cvm::error("Error: the CV " + name + 00033 " expects one or more nesting components.\n"); 00034 return; 00035 } else { 00036 // TODO: Maybe we can add an option to allow mixing scalar and vector types, 00037 // but that's a bit complicated so we just require a consistent type 00038 // of nesting CVs. 00039 x.type(cv[0]->value()); 00040 x.reset(); 00041 for (size_t i_cv = 1; i_cv < cv.size(); ++i_cv) { 00042 const auto type_i = cv[i_cv]->value().type(); 00043 if (type_i != x.type()) { 00044 cvm::error("Error: the type of sub-CVC " + cv[i_cv]->name + 00045 " is " + colvarvalue::type_desc(type_i) + ", which is " 00046 "different to the type of the first sub-CVC. Currently " 00047 "only sub-CVCs of the same type are supported to be " 00048 "nested.\n"); 00049 return; 00050 } 00051 } 00052 } 00053 use_explicit_gradients = true; 00054 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { 00055 if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { 00056 use_explicit_gradients = false; 00057 } 00058 } 00059 if (!use_explicit_gradients) { 00060 disable(f_cvc_explicit_gradient); 00061 } 00062 } 00063 00064 cvm::real colvar::linearCombination::getPolynomialFactorOfCVGradient(size_t i_cv) const { 00065 cvm::real factor_polynomial = 1.0; 00066 if (cv[i_cv]->value().type() == colvarvalue::type_scalar) { 00067 factor_polynomial = cv[i_cv]->sup_coeff * cv[i_cv]->sup_np * cvm::pow(cv[i_cv]->value().real_value, cv[i_cv]->sup_np - 1); 00068 } else { 00069 factor_polynomial = cv[i_cv]->sup_coeff; 00070 } 00071 return factor_polynomial; 00072 } 00073 00074 colvar::linearCombination::~linearCombination() { 00075 // Recall the steps we initialize the sub-CVCs: 00076 // 1. Lookup all sub-CVCs and then register the atom groups for sub-CVCs 00077 // in their constructors; 00078 // 2. Iterate over all sub-CVCs, get the pointers of their atom groups 00079 // groups, and register again in the parent (current) CVC. 00080 // That being said, the atom groups become children of the sub-CVCs at 00081 // first, and then become children of the parent CVC. 00082 // So, to destruct this class (parent CVC class), we need to remove the 00083 // dependencies of the atom groups to the parent CVC at first. 00084 remove_all_children(); 00085 // Then we remove the dependencies of the atom groups to the sub-CVCs 00086 // in their destructors. 00087 for (auto it = cv.begin(); it != cv.end(); ++it) { 00088 delete (*it); 00089 } 00090 // The last step is cleaning up the list of atom groups. 00091 atom_groups.clear(); 00092 } 00093 00094 void colvar::linearCombination::calc_value() { 00095 x.reset(); 00096 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { 00097 cv[i_cv]->calc_value(); 00098 colvarvalue current_cv_value(cv[i_cv]->value()); 00099 // polynomial combination allowed 00100 if (current_cv_value.type() == colvarvalue::type_scalar) { 00101 x += cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)); 00102 } else { 00103 x += cv[i_cv]->sup_coeff * current_cv_value; 00104 } 00105 } 00106 } 00107 00108 void colvar::linearCombination::calc_gradients() { 00109 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { 00110 cv[i_cv]->calc_gradients(); 00111 if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { 00112 cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); 00113 for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) { 00114 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { 00115 for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) { 00116 (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad; 00117 } 00118 } 00119 } 00120 } 00121 } 00122 } 00123 00124 void colvar::linearCombination::apply_force(colvarvalue const &force) { 00125 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { 00126 // If this CV us explicit gradients, then atomic gradients is already calculated 00127 // We can apply the force to atom groups directly 00128 if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { 00129 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { 00130 (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value); 00131 } 00132 } else { 00133 // Compute factors for polynomial combinations 00134 cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); 00135 colvarvalue cv_force = force.real_value * factor_polynomial; 00136 cv[i_cv]->apply_force(cv_force); 00137 } 00138 } 00139 } 00140 00141 colvar::customColvar::customColvar(std::string const &conf): linearCombination(conf) { 00142 use_custom_function = false; 00143 // code swipe from colvar::init_custom_function 00144 std::string expr_in, expr; 00145 size_t pos = 0; // current position in config string 00146 #ifdef LEPTON 00147 std::vector<Lepton::ParsedExpression> pexprs; 00148 Lepton::ParsedExpression pexpr; 00149 double *ref; 00150 #endif 00151 if (key_lookup(conf, "customFunction", &expr_in, &pos)) { 00152 #ifdef LEPTON 00153 use_custom_function = true; 00154 cvm::log("This colvar uses a custom function.\n"); 00155 do { 00156 expr = expr_in; 00157 if (cvm::debug()) 00158 cvm::log("Parsing expression \"" + expr + "\".\n"); 00159 try { 00160 pexpr = Lepton::Parser::parse(expr); 00161 pexprs.push_back(pexpr); 00162 } catch (...) { 00163 cvm::error("Error parsing expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR); 00164 } 00165 try { 00166 value_evaluators.push_back(new Lepton::CompiledExpression(pexpr.createCompiledExpression())); 00167 // Define variables for cvc values 00168 for (size_t i = 0; i < cv.size(); ++i) { 00169 for (size_t j = 0; j < cv[i]->value().size(); ++j) { 00170 std::string vn = cv[i]->name + (cv[i]->value().size() > 1 ? cvm::to_str(j+1) : ""); 00171 try { 00172 ref = &value_evaluators.back()->getVariableReference(vn); 00173 } catch (...) { 00174 ref = &dev_null; 00175 cvm::log("Warning: Variable " + vn + " is absent from expression \"" + expr + "\".\n"); 00176 } 00177 value_eval_var_refs.push_back(ref); 00178 } 00179 } 00180 } catch (...) { 00181 cvm::error("Error compiling expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR); 00182 } 00183 } while (key_lookup(conf, "customFunction", &expr_in, &pos)); 00184 // Now define derivative with respect to each scalar sub-component 00185 for (size_t i = 0; i < cv.size(); ++i) { 00186 for (size_t j = 0; j < cv[i]->value().size(); ++j) { 00187 std::string vn = cv[i]->name + (cv[i]->value().size() > 1 ? cvm::to_str(j+1) : ""); 00188 for (size_t c = 0; c < pexprs.size(); ++c) { 00189 gradient_evaluators.push_back(new Lepton::CompiledExpression(pexprs[c].differentiate(vn).createCompiledExpression())); 00190 for (size_t k = 0; k < cv.size(); ++k) { 00191 for (size_t l = 0; l < cv[k]->value().size(); l++) { 00192 std::string vvn = cv[k]->name + (cv[k]->value().size() > 1 ? cvm::to_str(l+1) : ""); 00193 try { 00194 ref = &gradient_evaluators.back()->getVariableReference(vvn); 00195 } catch (...) { 00196 cvm::log("Warning: Variable " + vvn + " is absent from derivative of \"" + expr + "\" wrt " + vn + ".\n"); 00197 ref = &dev_null; 00198 } 00199 grad_eval_var_refs.push_back(ref); 00200 } 00201 } 00202 } 00203 } 00204 } 00205 if (value_evaluators.size() == 0) { 00206 cvm::error("Error: no custom function defined.\n", COLVARS_INPUT_ERROR); 00207 } 00208 if (value_evaluators.size() != 1) { 00209 x.type(colvarvalue::type_vector); 00210 } else { 00211 x.type(colvarvalue::type_scalar); 00212 } 00213 #else 00214 cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" 00215 "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", 00216 COLVARS_INPUT_ERROR); 00217 #endif 00218 } else { 00219 cvm::log("Warning: no customFunction specified.\n"); 00220 cvm::log("Warning: use linear combination instead.\n"); 00221 } 00222 } 00223 00224 colvar::customColvar::~customColvar() { 00225 #ifdef LEPTON 00226 for (size_t i = 0; i < value_evaluators.size(); ++i) { 00227 if (value_evaluators[i] != nullptr) delete value_evaluators[i]; 00228 } 00229 for (size_t i = 0; i < gradient_evaluators.size(); ++i) { 00230 if (gradient_evaluators[i] != nullptr) delete gradient_evaluators[i]; 00231 } 00232 #endif 00233 } 00234 00235 void colvar::customColvar::calc_value() { 00236 if (!use_custom_function) { 00237 colvar::linearCombination::calc_value(); 00238 } else { 00239 #ifdef LEPTON 00240 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { 00241 cv[i_cv]->calc_value(); 00242 } 00243 x.reset(); 00244 size_t l = 0; 00245 for (size_t i = 0; i < x.size(); ++i) { 00246 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { 00247 const colvarvalue& current_cv_value = cv[i_cv]->value(); 00248 for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { 00249 if (current_cv_value.type() == colvarvalue::type_scalar) { 00250 *(value_eval_var_refs[l++]) = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)); 00251 } else { 00252 *(value_eval_var_refs[l++]) = cv[i_cv]->sup_coeff * current_cv_value[j_elem]; 00253 } 00254 } 00255 } 00256 x[i] = value_evaluators[i]->evaluate(); 00257 } 00258 #else 00259 cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" 00260 "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", 00261 COLVARS_INPUT_ERROR); 00262 #endif 00263 } 00264 } 00265 00266 void colvar::customColvar::calc_gradients() { 00267 if (!use_custom_function) { 00268 colvar::linearCombination::calc_gradients(); 00269 } else { 00270 #ifdef LEPTON 00271 size_t r = 0; // index in the vector of variable references 00272 size_t e = 0; // index of the gradient evaluator 00273 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { // for each CV 00274 cv[i_cv]->calc_gradients(); 00275 if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { 00276 const colvarvalue& current_cv_value = cv[i_cv]->value(); 00277 const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); 00278 for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { // for each element in this CV 00279 for (size_t c = 0; c < x.size(); ++c) { // for each custom function expression 00280 for (size_t k = 0; k < cv.size(); ++k) { // this is required since we need to feed all CV values to this expression 00281 const cvm::real factor_polynomial_k = getPolynomialFactorOfCVGradient(k); 00282 for (size_t l = 0; l < cv[k]->value().size(); ++l) { 00283 *(grad_eval_var_refs[r++]) = factor_polynomial_k * cv[k]->value()[l]; 00284 } 00285 } 00286 const double expr_grad = gradient_evaluators[e++]->evaluate(); 00287 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { 00288 for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) { 00289 (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = expr_grad * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad; 00290 } 00291 } 00292 } 00293 } 00294 } 00295 } 00296 #else 00297 cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" 00298 "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", 00299 COLVARS_INPUT_ERROR); 00300 #endif 00301 } 00302 } 00303 00304 void colvar::customColvar::apply_force(colvarvalue const &force) { 00305 if (!use_custom_function) { 00306 colvar::linearCombination::apply_force(force); 00307 } else { 00308 #ifdef LEPTON 00309 size_t r = 0; // index in the vector of variable references 00310 size_t e = 0; // index of the gradient evaluator 00311 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { 00312 // If this CV us explicit gradients, then atomic gradients is already calculated 00313 // We can apply the force to atom groups directly 00314 if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { 00315 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { 00316 (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value); 00317 } 00318 } else { 00319 const colvarvalue& current_cv_value = cv[i_cv]->value(); 00320 colvarvalue cv_force(current_cv_value.type()); 00321 const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); 00322 for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { 00323 for (size_t c = 0; c < x.size(); ++c) { 00324 for (size_t k = 0; k < cv.size(); ++k) { 00325 const cvm::real factor_polynomial_k = getPolynomialFactorOfCVGradient(k); 00326 for (size_t l = 0; l < cv[k]->value().size(); ++l) { 00327 *(grad_eval_var_refs[r++]) = factor_polynomial_k * cv[k]->value()[l]; 00328 } 00329 } 00330 cv_force[j_elem] += factor_polynomial * gradient_evaluators[e++]->evaluate() * force.real_value; 00331 } 00332 } 00333 cv[i_cv]->apply_force(cv_force); 00334 } 00335 } 00336 #else 00337 cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" 00338 "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", 00339 COLVARS_INPUT_ERROR); 00340 #endif 00341 } 00342 } 00343 00344 #endif // __cplusplus >= 201103L