00001 #if (__cplusplus >= 201103L)
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "colvarcomp.h"
00011
00012 colvar::linearCombination::linearCombination(std::string const &conf): cvc(conf) {
00013
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
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
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
00037
00038
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
00076
00077
00078
00079
00080
00081
00082
00083
00084 remove_all_children();
00085
00086
00087 for (auto it = cv.begin(); it != cv.end(); ++it) {
00088 delete (*it);
00089 }
00090
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
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
00127
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
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
00144 std::string expr_in, expr;
00145 size_t pos = 0;
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
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
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;
00272 size_t e = 0;
00273 for (size_t i_cv = 0; i_cv < cv.size(); ++i_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) {
00279 for (size_t c = 0; c < x.size(); ++c) {
00280 for (size_t k = 0; k < cv.size(); ++k) {
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;
00310 size_t e = 0;
00311 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00312
00313
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