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 #ifndef COLVAR_H 00011 #define COLVAR_H 00012 00013 #include <iostream> 00014 00015 #if (__cplusplus >= 201103L) 00016 #include <map> 00017 #include <functional> 00018 #endif 00019 00020 #include "colvarmodule.h" 00021 #include "colvarvalue.h" 00022 #include "colvarparse.h" 00023 #include "colvardeps.h" 00024 00025 #ifdef LEPTON 00026 #include "Lepton.h" // for runtime custom expressions 00027 #endif 00028 00053 00054 class colvar : public colvarparse, public colvardeps { 00055 00056 public: 00057 00059 std::string name; 00060 00062 colvarvalue const & value() const; 00063 00065 colvarvalue const & actual_value() const; 00066 00068 colvarvalue const & run_ave() const; 00069 00071 cvm::real const & force_constant() const; 00072 00074 colvarvalue const & velocity() const; 00075 00084 colvarvalue const & total_force() const; 00085 00094 cvm::real width; 00095 00097 static std::vector<feature *> cv_features; 00098 00100 virtual const std::vector<feature *> &features() const 00101 { 00102 return cv_features; 00103 } 00104 virtual std::vector<feature *> &modify_features() 00105 { 00106 return cv_features; 00107 } 00108 static void delete_features() { 00109 for (size_t i=0; i < cv_features.size(); i++) { 00110 delete cv_features[i]; 00111 } 00112 cv_features.clear(); 00113 } 00114 00118 void do_feature_side_effects(int id); 00119 00121 std::vector<colvarbias *> biases; 00122 00123 protected: 00124 00125 00126 /* 00127 extended: 00128 H = H_{0} + \sum_{i} 1/2*\lambda*(S_i(x(t))-s_i(t))^2 \\ 00129 + \sum_{i} 1/2*m_i*(ds_i(t)/dt)^2 \\ 00130 + \sum_{t'<t} W * exp(-1/2*\sum_{i} (s_i(t')-s_i(t))^2/(\delta{}s_i)^2) \\ 00131 + \sum_{w} (\sum_{i}c_{w,i}s_i(t) - D_w)^M/(\Sigma_w)^M 00132 00133 normal: 00134 H = H_{0} + \sum_{t'<t} W * exp(-1/2*\sum_{i} (S_i(x(t'))-S_i(x(t)))^2/(\delta{}S_i)^2) \\ 00135 + \sum_{w} (\sum_{i}c_{w,i}S_i(t) - D_w)^M/(\Sigma_w)^M 00136 00137 output: 00138 H = H_{0} (only output S(x), no forces) 00139 00140 Here: 00141 S(x(t)) = x 00142 s(t) = x_ext 00143 DS = Ds = delta 00144 */ 00145 00146 00148 colvarvalue x; 00149 00150 // TODO: implement functionality to treat these 00151 // /// Vector of individual values from CVCs 00152 // colvarvalue x_cvc; 00153 00154 // /// Jacobian matrix of individual values from CVCs 00155 // colvarvalue dx_cvc; 00156 00158 colvarvalue x_reported; 00159 00161 colvarvalue v_fdiff; 00162 00163 inline colvarvalue fdiff_velocity(colvarvalue const &xold, 00164 colvarvalue const &xnew) 00165 { 00166 // using the gradient of the square distance to calculate the 00167 // velocity (non-scalar variables automatically taken into 00168 // account) 00169 cvm::real dt = cvm::dt(); 00170 return ( ( (dt > 0.0) ? (1.0/dt) : 1.0 ) * 00171 0.5 * dist2_lgrad(xnew, xold) ); 00172 } 00173 00175 colvarvalue v_reported; 00176 00177 // Options for extended_lagrangian 00179 colvarvalue x_ext; 00181 colvarvalue prev_x_ext; 00183 colvarvalue v_ext; 00185 colvarvalue prev_v_ext; 00187 cvm::real ext_mass; 00189 cvm::real ext_force_k; 00191 cvm::real ext_gamma; 00193 cvm::real ext_sigma; 00194 00196 colvarvalue fr; 00197 00199 colvarvalue fj; 00200 00202 colvarvalue ft_reported; 00203 00204 public: 00205 00206 00209 colvarvalue fb; 00210 00212 colvarvalue fb_actual; 00213 00217 colvarvalue f; 00218 00220 colvarvalue f_old; 00221 00224 colvarvalue ft; 00225 00227 cvm::real period; 00228 00230 cvm::real wrap_center; 00231 00234 bool expand_boundaries; 00235 00237 colvarvalue lower_boundary; 00238 00240 colvarvalue upper_boundary; 00241 00243 bool periodic_boundaries() const; 00244 00246 bool periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const; 00247 00248 00250 colvar(); 00251 00253 int init(std::string const &conf); 00254 00256 int init_components(std::string const &conf); 00257 00259 int init_custom_function(std::string const &conf); 00260 00262 int init_grid_parameters(std::string const &conf); 00263 00265 int init_extended_Lagrangian(std::string const &conf); 00266 00268 int init_output_flags(std::string const &conf); 00269 00271 virtual int init_dependencies(); 00272 00273 private: 00275 template<typename def_class_name> int init_components_type(std::string const & conf, 00276 char const *def_desc, 00277 char const *def_config_key); 00278 #if (__cplusplus >= 201103L) 00282 int init_components_type_from_global_map(const std::string& conf, 00283 const char* def_config_key); 00284 #endif 00285 00286 public: 00287 00289 void setup(); 00290 00292 ~colvar(); 00293 00294 00296 int calc(); 00297 00299 int end_of_step(); 00300 00304 int calc_cvcs(int first = 0, size_t num_cvcs = 0); 00305 00307 int check_cvc_range(int first_cvc, size_t num_cvcs); 00308 00310 int calc_cvc_values(int first, size_t num_cvcs); 00312 int calc_cvc_gradients(int first, size_t num_cvcs); 00314 int calc_cvc_total_force(int first, size_t num_cvcs); 00316 int calc_cvc_Jacobians(int first, size_t num_cvcs); 00317 00319 int collect_cvc_data(); 00320 00322 int collect_cvc_values(); 00324 int collect_cvc_gradients(); 00326 int collect_cvc_total_forces(); 00328 int collect_cvc_Jacobians(); 00330 int calc_colvar_properties(); 00331 00333 inline colvarvalue const applied_force() const 00334 { 00335 if (is_enabled(f_cv_extended_Lagrangian)) { 00336 return fr; 00337 } 00338 return f; 00339 } 00340 00342 void reset_bias_force(); 00343 00345 void add_bias_force(colvarvalue const &force); 00346 00348 void add_bias_force_actual_value(colvarvalue const &force); 00349 00354 cvm::real update_forces_energy(); 00355 00357 void update_extended_Lagrangian(); 00358 00361 void communicate_forces(); 00362 00364 int set_cvc_flags(std::vector<bool> const &flags); 00365 00367 int update_cvc_flags(); 00368 00370 int update_cvc_config(std::vector<std::string> const &confs); 00371 00373 int cvc_param_exists(std::string const ¶m_name); 00374 00376 cvm::real get_cvc_param(std::string const ¶m_name); 00377 00379 void const *get_cvc_param_ptr(std::string const ¶m_name); 00380 00382 colvarvalue const *get_cvc_param_grad(std::string const ¶m_name); 00383 00385 int set_cvc_param(std::string const ¶m_name, void const *new_value); 00386 00387 protected: 00388 00390 size_t n_active_cvcs; 00391 00393 cvm::real active_cvc_square_norm; 00394 00396 void update_active_cvc_square_norm(); 00397 00399 cvm::step_number prev_timestep; 00400 00401 public: 00402 00404 inline size_t num_dimensions() const 00405 { 00406 return value().size(); 00407 } 00408 00410 inline size_t num_cvcs() const 00411 { 00412 return cvcs.size(); 00413 } 00414 00417 inline size_t num_active_cvcs() const 00418 { 00419 return n_active_cvcs; 00420 } 00421 00426 cvm::real dist2(colvarvalue const &x1, 00427 colvarvalue const &x2) const; 00428 00433 colvarvalue dist2_lgrad(colvarvalue const &x1, 00434 colvarvalue const &x2) const; 00435 00440 colvarvalue dist2_rgrad(colvarvalue const &x1, 00441 colvarvalue const &x2) const; 00442 00447 void wrap(colvarvalue &x_unwrapped) const; 00448 00450 int parse_analysis(std::string const &conf); 00451 00453 int analyze(); 00454 00456 std::istream & read_traj(std::istream &is); 00457 00459 std::ostream & write_traj(std::ostream &os); 00461 std::ostream & write_traj_label(std::ostream &os); 00462 00464 std::istream & read_state(std::istream &is); 00466 std::ostream & write_state(std::ostream &os); 00467 00469 int write_output_files(); 00470 00471 protected: 00472 00474 colvarvalue x_old; 00475 00477 colvarvalue x_restart; 00478 00480 bool after_restart; 00481 00484 std::list< std::list<colvarvalue> > acf_x_history, acf_v_history; 00487 std::list< std::list<colvarvalue> >::iterator acf_x_history_p, acf_v_history_p; 00488 00490 std::list< std::list<colvarvalue> > x_history; 00493 std::list< std::list<colvarvalue> >::iterator x_history_p; 00494 00497 std::string acf_colvar_name; 00499 size_t acf_length; 00501 size_t acf_offset; 00503 size_t acf_stride; 00505 size_t acf_nframes; 00507 bool acf_normalize; 00509 std::vector<cvm::real> acf; 00511 std::string acf_outfile; 00512 00514 enum acf_type_e { 00516 acf_notset, 00518 acf_vel, 00520 acf_coor, 00523 acf_p2coor 00524 }; 00525 00527 acf_type_e acf_type; 00528 00530 void calc_vel_acf(std::list<colvarvalue> &v_history, 00531 colvarvalue const &v); 00532 00535 void calc_coor_acf(std::list<colvarvalue> &x_history, 00536 colvarvalue const &x); 00537 00540 void calc_p2coor_acf(std::list<colvarvalue> &x_history, 00541 colvarvalue const &x); 00542 00544 int calc_acf(); 00546 int write_acf(std::ostream &os); 00547 00549 size_t runave_length; 00551 size_t runave_stride; 00553 std::string runave_outfile; 00555 colvarvalue runave; 00557 cvm::real runave_variance; 00558 00560 int calc_runave(); 00561 00563 cvm::real kinetic_energy; 00565 cvm::real potential_energy; 00566 00567 public: 00568 00569 // collective variable component base class 00570 class cvc; 00571 00572 // list of available collective variable components 00573 00574 // scalar colvar components 00575 class distance; 00576 class distance_z; 00577 class distance_xy; 00578 class polar_theta; 00579 class polar_phi; 00580 class distance_inv; 00581 class distance_pairs; 00582 class dipole_magnitude; 00583 class angle; 00584 class dipole_angle; 00585 class dihedral; 00586 class coordnum; 00587 class selfcoordnum; 00588 class groupcoordnum; 00589 class h_bond; 00590 class rmsd; 00591 class orientation_angle; 00592 class orientation_proj; 00593 class tilt; 00594 class spin_angle; 00595 class gyration; 00596 class inertia; 00597 class inertia_z; 00598 class eigenvector; 00599 class alpha_dihedrals; 00600 class alpha_angles; 00601 class dihedPC; 00602 class alch_lambda; 00603 class alch_Flambda; 00604 class componentDisabled; 00605 class CartesianBasedPath; 00606 class gspath; 00607 class gzpath; 00608 class linearCombination; 00609 class CVBasedPath; 00610 class gspathCV; 00611 class gzpathCV; 00612 class aspathCV; 00613 class azpathCV; 00614 class euler_phi; 00615 class euler_psi; 00616 class euler_theta; 00617 class neuralNetwork; 00618 class customColvar; 00619 00620 // non-scalar components 00621 class distance_vec; 00622 class distance_dir; 00623 class cartesian; 00624 class orientation; 00625 00626 // components that do not handle any atoms directly 00627 class map_total; 00628 00630 #if (__cplusplus >= 201103L) 00631 00632 static const std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>>& get_global_cvc_map() { 00633 return global_cvc_map; 00634 } 00635 #endif 00636 00638 static bool compare_cvc(const colvar::cvc* const i, const colvar::cvc* const j); 00639 00640 protected: 00641 00643 std::vector<cvc *> cvcs; 00644 00646 std::vector<bool> cvc_flags; 00647 00650 void build_atom_list(void); 00651 00653 std::string scripted_function; 00654 00657 std::vector<const colvarvalue *> sorted_cvc_values; 00658 00659 #ifdef LEPTON 00660 00661 std::vector<Lepton::CompiledExpression *> value_evaluators; 00662 00664 std::vector<Lepton::CompiledExpression *> gradient_evaluators; 00665 00667 std::vector<double *> value_eval_var_refs; 00668 std::vector<double *> grad_eval_var_refs; 00669 00671 double dev_null; 00672 #endif 00673 00674 #if (__cplusplus >= 201103L) 00675 00676 static std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>> global_cvc_map; 00677 #endif 00678 00680 std::vector<int> volmap_ids_; 00681 00682 public: 00683 00685 std::vector<int> atom_ids; 00686 00690 std::vector<cvm::rvector> atomic_gradients; 00691 00693 virtual std::vector<std::vector<int> > get_atom_lists(); 00694 00696 std::vector<int> const &get_volmap_ids(); 00697 00698 }; 00699 00700 00701 inline cvm::real const & colvar::force_constant() const 00702 { 00703 return ext_force_k; 00704 } 00705 00706 00707 inline colvarvalue const & colvar::value() const 00708 { 00709 return x_reported; 00710 } 00711 00712 00713 inline colvarvalue const & colvar::actual_value() const 00714 { 00715 return x; 00716 } 00717 00718 00719 inline colvarvalue const & colvar::run_ave() const 00720 { 00721 return runave; 00722 } 00723 00724 00725 inline colvarvalue const & colvar::velocity() const 00726 { 00727 return v_reported; 00728 } 00729 00730 00731 inline colvarvalue const & colvar::total_force() const 00732 { 00733 return ft_reported; 00734 } 00735 00736 00737 inline void colvar::add_bias_force(colvarvalue const &force) 00738 { 00739 check_enabled(f_cv_gradient, 00740 std::string("applying a force to the variable \""+name+"\"")); 00741 if (cvm::debug()) { 00742 cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n"); 00743 } 00744 fb += force; 00745 } 00746 00747 00748 inline void colvar::add_bias_force_actual_value(colvarvalue const &force) 00749 { 00750 if (cvm::debug()) { 00751 cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n"); 00752 } 00753 fb_actual += force; 00754 } 00755 00756 00757 inline void colvar::reset_bias_force() { 00758 fb.type(value()); 00759 fb.reset(); 00760 fb_actual.type(value()); 00761 fb_actual.reset(); 00762 } 00763 00764 #endif 00765