00001 // -*- c++ -*- 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 <iostream> 00011 #include <iomanip> 00012 #include <cstdlib> 00013 00014 #include "colvarmodule.h" 00015 #include "colvarproxy.h" 00016 #include "colvarbias.h" 00017 #include "colvarbias_alb.h" 00018 00019 #ifdef _MSC_VER 00020 #if _MSC_VER <= 1700 00021 #define copysign(A,B) _copysign(A,B) 00022 double fmax(double A, double B) { return ( A > B ? A : B ); } 00023 double fmin(double A, double B) { return ( A < B ? A : B ); } 00024 #endif 00025 #endif 00026 00027 /* Note about nomenclature. Force constant is called a coupling 00028 * constant here to emphasize its changing in the code. Outwards, 00029 * everything is called a force constant to keep it consistent with 00030 * the rest of colvars. 00031 * 00032 */ 00033 00034 colvarbias_alb::colvarbias_alb(char const *key) 00035 : colvarbias(key), update_calls(0), b_equilibration(true) 00036 { 00037 } 00038 00039 00040 int colvarbias_alb::init(std::string const &conf) 00041 { 00042 colvarproxy *proxy = cvm::main()->proxy; 00043 colvarbias::init(conf); 00044 cvm::main()->cite_feature("ALB colvar bias implementation"); 00045 00046 enable(f_cvb_scalar_variables); 00047 00048 size_t i; 00049 00050 // get the initial restraint centers 00051 colvar_centers.resize(num_variables()); 00052 00053 means.resize(num_variables()); 00054 ssd.resize(num_variables()); //sum of squares of differences from mean 00055 00056 //setup force vectors 00057 max_coupling_range.resize(num_variables()); 00058 max_coupling_rate.resize(num_variables()); 00059 coupling_accum.resize(num_variables()); 00060 set_coupling.resize(num_variables()); 00061 current_coupling.resize(num_variables()); 00062 coupling_rate.resize(num_variables()); 00063 00064 enable(f_cvb_apply_force); 00065 00066 for (i = 0; i < num_variables(); i++) { 00067 colvar_centers[i].type(colvars[i]->value()); 00068 //zero moments 00069 means[i] = ssd[i] = 0; 00070 00071 //zero force some of the force vectors that aren't initialized 00072 coupling_accum[i] = current_coupling[i] = 0; 00073 00074 } 00075 if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { 00076 for (i = 0; i < num_variables(); i++) { 00077 colvar_centers[i].apply_constraints(); 00078 } 00079 } else { 00080 colvar_centers.clear(); 00081 cvm::error("Error: must define the initial centers of adaptive linear bias .\n"); 00082 } 00083 00084 if (colvar_centers.size() != num_variables()) 00085 cvm::error("Error: number of centers does not match " 00086 "that of collective variables.\n"); 00087 00088 if (!get_keyval(conf, "UpdateFrequency", update_freq, 0)) 00089 cvm::error("Error: must set updateFrequency for adaptive linear bias.\n"); 00090 00091 //we split the time between updating and equilibrating 00092 update_freq /= 2; 00093 00094 if (update_freq <= 1) 00095 cvm::error("Error: must set updateFrequency to greater than 2.\n"); 00096 00097 enable(f_cvb_history_dependent); 00098 00099 get_keyval(conf, "outputCenters", b_output_centers, false); 00100 get_keyval(conf, "outputGradient", b_output_grad, false); 00101 get_keyval(conf, "outputCoupling", b_output_coupling, true); 00102 get_keyval(conf, "hardForceRange", b_hard_coupling_range, true); 00103 00104 //initial guess 00105 if (!get_keyval(conf, "forceConstant", set_coupling, set_coupling)) 00106 for (i =0 ; i < num_variables(); i++) 00107 set_coupling[i] = 0.; 00108 00109 //how we're going to increase to that point 00110 for (i = 0; i < num_variables(); i++) 00111 coupling_rate[i] = (set_coupling[i] - current_coupling[i]) / update_freq; 00112 00113 00114 if (!get_keyval(conf, "forceRange", max_coupling_range, max_coupling_range)) { 00115 //set to default 00116 for (i = 0; i < num_variables(); i++) { 00117 if (proxy->target_temperature() > 0.0) { 00118 max_coupling_range[i] = 3 * proxy->target_temperature() * 00119 proxy->boltzmann(); 00120 } else { 00121 max_coupling_range[i] = 3 * proxy->boltzmann(); 00122 } 00123 } 00124 } 00125 00126 if (!get_keyval(conf, "rateMax", max_coupling_rate, max_coupling_rate)) { 00127 //set to default 00128 for (i = 0; i < num_variables(); i++) { 00129 max_coupling_rate[i] = max_coupling_range[i] / (10 * update_freq); 00130 } 00131 } 00132 00133 00134 if (cvm::debug()) 00135 cvm::log(" bias.\n"); 00136 00137 return COLVARS_OK; 00138 } 00139 00140 00141 colvarbias_alb::~colvarbias_alb() 00142 { 00143 } 00144 00145 00146 int colvarbias_alb::update() 00147 { 00148 colvarproxy *proxy = cvm::main()->proxy; 00149 00150 bias_energy = 0.0; 00151 update_calls++; 00152 00153 if (cvm::debug()) 00154 cvm::log("Updating the adaptive linear bias \""+this->name+"\".\n"); 00155 00156 //log the moments of the CVs 00157 // Force and energy calculation 00158 bool finished_equil_flag = 1; 00159 cvm::real delta; 00160 for (size_t i = 0; i < num_variables(); i++) { 00161 colvar_forces[i] = -1.0 * restraint_force(restraint_convert_k(current_coupling[i], colvars[i]->width), 00162 colvars[i], 00163 colvar_centers[i]); 00164 bias_energy += restraint_potential(restraint_convert_k(current_coupling[i], colvars[i]->width), 00165 colvars[i], 00166 colvar_centers[i]); 00167 00168 if (!b_equilibration) { 00169 //Welford, West, and Hanso online variance method 00170 00171 delta = static_cast<cvm::real>(colvars[i]->value()) - means[i]; 00172 means[i] += delta / update_calls; 00173 ssd[i] += delta * (static_cast<cvm::real>(colvars[i]->value()) - means[i]); 00174 00175 } else { 00176 //check if we've reached the setpoint 00177 cvm::real const coupling_diff = current_coupling[i] - set_coupling[i]; 00178 if ((coupling_rate[i] == 0) || 00179 ((coupling_diff*coupling_diff) < (coupling_rate[i]*coupling_rate[i]))) { 00180 finished_equil_flag &= 1; //we continue equilibrating as long as we haven't reached all the set points 00181 } 00182 else { 00183 current_coupling[i] += coupling_rate[i]; 00184 finished_equil_flag = 0; 00185 } 00186 00187 00188 //update max_coupling_range 00189 if (!b_hard_coupling_range && fabs(current_coupling[i]) > max_coupling_range[i]) { 00190 std::ostringstream logStream; 00191 logStream << "Coupling constant for " 00192 << colvars[i]->name 00193 << " has exceeded coupling range of " 00194 << max_coupling_range[i] 00195 << ".\n"; 00196 00197 max_coupling_range[i] *= 1.25; 00198 logStream << "Expanding coupling range to " << max_coupling_range[i] << ".\n"; 00199 cvm::log(logStream.str()); 00200 } 00201 00202 00203 } 00204 } 00205 00206 if (b_equilibration && finished_equil_flag) { 00207 b_equilibration = false; 00208 update_calls = 0; 00209 } 00210 00211 00212 //now we update coupling constant, if necessary 00213 if (!b_equilibration && update_calls == update_freq) { 00214 00215 //use estimated variance to take a step 00216 cvm::real step_size = 0; 00217 cvm::real temp; 00218 00219 //reset means and sum of squares of differences 00220 for (size_t i = 0; i < num_variables(); i++) { 00221 00222 temp = 2. * (means[i] / (static_cast<cvm::real> (colvar_centers[i])) - 1) * ssd[i] / (update_calls - 1); 00223 00224 if (proxy->target_temperature() > 0.0) { 00225 step_size = temp / (proxy->target_temperature() * proxy->boltzmann()); 00226 } else { 00227 step_size = temp / proxy->boltzmann(); 00228 } 00229 00230 means[i] = 0; 00231 ssd[i] = 0; 00232 00233 //stochastic if we do that update or not 00234 if (num_variables() == 1 || rand() < RAND_MAX / ((int) num_variables())) { 00235 coupling_accum[i] += step_size * step_size; 00236 current_coupling[i] = set_coupling[i]; 00237 set_coupling[i] += max_coupling_range[i] / sqrt(coupling_accum[i]) * step_size; 00238 coupling_rate[i] = (set_coupling[i] - current_coupling[i]) / update_freq; 00239 //set to the minimum rate and then put the sign back on it 00240 coupling_rate[i] = copysign(fmin(fabs(coupling_rate[i]), max_coupling_rate[i]), coupling_rate[i]); 00241 } else { 00242 coupling_rate[i] = 0; 00243 } 00244 00245 } 00246 00247 update_calls = 0; 00248 b_equilibration = true; 00249 00250 } 00251 00252 return COLVARS_OK; 00253 } 00254 00255 00256 int colvarbias_alb::set_state_params(std::string const &conf) 00257 { 00258 int error_code = colvarbias::set_state_params(conf); 00259 00260 if (error_code != COLVARS_OK) { 00261 return error_code; 00262 } 00263 00264 if (!get_keyval(conf, "setCoupling", set_coupling)) 00265 cvm::error("Error: current setCoupling is missing from the restart.\n"); 00266 00267 if (!get_keyval(conf, "currentCoupling", current_coupling)) 00268 cvm::error("Error: current setCoupling is missing from the restart.\n"); 00269 00270 if (!get_keyval(conf, "maxCouplingRange", max_coupling_range)) 00271 cvm::error("Error: maxCouplingRange is missing from the restart.\n"); 00272 00273 if (!get_keyval(conf, "couplingRate", coupling_rate)) 00274 cvm::error("Error: current setCoupling is missing from the restart.\n"); 00275 00276 if (!get_keyval(conf, "couplingAccum", coupling_accum)) 00277 cvm::error("Error: couplingAccum is missing from the restart.\n"); 00278 00279 if (!get_keyval(conf, "mean", means)) 00280 cvm::error("Error: current mean is missing from the restart.\n"); 00281 00282 if (!get_keyval(conf, "ssd", ssd)) 00283 cvm::error("Error: current ssd is missing from the restart.\n"); 00284 00285 if (!get_keyval(conf, "updateCalls", update_calls)) 00286 cvm::error("Error: current updateCalls is missing from the restart.\n"); 00287 00288 if (!get_keyval(conf, "b_equilibration", b_equilibration)) 00289 cvm::error("Error: current updateCalls is missing from the restart.\n"); 00290 00291 return COLVARS_OK; 00292 } 00293 00294 00295 std::string const colvarbias_alb::get_state_params() const 00296 { 00297 std::ostringstream os; 00298 os << " setCoupling "; 00299 size_t i; 00300 for (i = 0; i < num_variables(); i++) { 00301 os << std::setprecision(cvm::en_prec) 00302 << std::setw(cvm::en_width) << set_coupling[i] << "\n"; 00303 } 00304 os << " currentCoupling "; 00305 for (i = 0; i < num_variables(); i++) { 00306 os << std::setprecision(cvm::en_prec) 00307 << std::setw(cvm::en_width) << current_coupling[i] << "\n"; 00308 } 00309 os << " maxCouplingRange "; 00310 for (i = 0; i < num_variables(); i++) { 00311 os << std::setprecision(cvm::en_prec) 00312 << std::setw(cvm::en_width) << max_coupling_range[i] << "\n"; 00313 } 00314 os << " couplingRate "; 00315 for (i = 0; i < num_variables(); i++) { 00316 os << std::setprecision(cvm::en_prec) 00317 << std::setw(cvm::en_width) << coupling_rate[i] << "\n"; 00318 } 00319 os << " couplingAccum "; 00320 for (i = 0; i < num_variables(); i++) { 00321 os << std::setprecision(cvm::en_prec) 00322 << std::setw(cvm::en_width) << coupling_accum[i] << "\n"; 00323 } 00324 os << " mean "; 00325 for (i = 0; i < num_variables(); i++) { 00326 os << std::setprecision(cvm::en_prec) 00327 << std::setw(cvm::en_width) << means[i] << "\n"; 00328 } 00329 os << " ssd "; 00330 for (i = 0; i < num_variables(); i++) { 00331 os << std::setprecision(cvm::en_prec) 00332 << std::setw(cvm::en_width) << ssd[i] << "\n"; 00333 } 00334 os << " updateCalls " << update_calls << "\n"; 00335 if (b_equilibration) 00336 os << " b_equilibration yes\n"; 00337 else 00338 os << " b_equilibration no\n"; 00339 00340 return os.str(); 00341 } 00342 00343 00344 std::ostream & colvarbias_alb::write_traj_label(std::ostream &os) 00345 { 00346 os << " "; 00347 00348 if (b_output_energy) 00349 os << " E_" 00350 << cvm::wrap_string(this->name, cvm::en_width-2); 00351 00352 if (b_output_coupling) 00353 for (size_t i = 0; i < current_coupling.size(); i++) { 00354 os << " ForceConst_" << i 00355 <<std::setw(cvm::en_width - 6 - (i / 10 + 1)) 00356 << ""; 00357 } 00358 00359 if (b_output_grad) 00360 for (size_t i = 0; i < means.size(); i++) { 00361 os << "Grad_" 00362 << cvm::wrap_string(colvars[i]->name, cvm::cv_width - 4); 00363 } 00364 00365 if (b_output_centers) 00366 for (size_t i = 0; i < num_variables(); i++) { 00367 size_t const this_cv_width = (colvars[i]->value()).output_width(cvm::cv_width); 00368 os << " x0_" 00369 << cvm::wrap_string(colvars[i]->name, this_cv_width-3); 00370 } 00371 00372 return os; 00373 } 00374 00375 00376 std::ostream & colvarbias_alb::write_traj(std::ostream &os) 00377 { 00378 os << " "; 00379 00380 if (b_output_energy) 00381 os << " " 00382 << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width) 00383 << bias_energy; 00384 00385 if (b_output_coupling) 00386 for (size_t i = 0; i < current_coupling.size(); i++) { 00387 os << " " 00388 << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width) 00389 << current_coupling[i]; 00390 } 00391 00392 00393 if (b_output_centers) 00394 for (size_t i = 0; i < num_variables(); i++) { 00395 os << " " 00396 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) 00397 << colvar_centers[i]; 00398 } 00399 00400 if (b_output_grad) 00401 for (size_t i = 0; i < means.size(); i++) { 00402 os << " " 00403 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) 00404 << -2.0 * (means[i] / (static_cast<cvm::real>(colvar_centers[i])) - 1) * ssd[i] / (fmax(update_calls, 2.0) - 1); 00405 00406 } 00407 00408 return os; 00409 } 00410 00411 00412 cvm::real colvarbias_alb::restraint_potential(cvm::real k, 00413 colvar const *x, 00414 colvarvalue const &xcenter) const 00415 { 00416 return k * (x->value() - xcenter); 00417 } 00418 00419 00420 colvarvalue colvarbias_alb::restraint_force(cvm::real k, 00421 colvar const * /* x */, 00422 colvarvalue const & /* xcenter */) const 00423 { 00424 return k; 00425 } 00426 00427 00428 cvm::real colvarbias_alb::restraint_convert_k(cvm::real k, 00429 cvm::real dist_measure) const 00430 { 00431 return k / dist_measure; 00432 } 00433