Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

colvargrid.h

Go to the documentation of this file.
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 COLVARGRID_H
00011 #define COLVARGRID_H
00012 
00013 #include <iostream>
00014 #include <iomanip>
00015 
00016 #include "colvar.h"
00017 #include "colvarmodule.h"
00018 #include "colvarvalue.h"
00019 #include "colvarparse.h"
00020 
00025 template <class T> class colvar_grid : public colvarparse {
00026 
00027 protected:
00028 
00030   size_t nd;
00031 
00033   std::vector<int> nx;
00034 
00036   std::vector<int> nxc;
00037 
00040   size_t mult;
00041 
00043   size_t nt;
00044 
00046   std::vector<T> data;
00047 
00049   std::vector<size_t> new_data;
00050 
00052   std::vector<colvar *> cv;
00053 
00055   std::vector<bool> use_actual_value;
00056 
00058   inline size_t address(std::vector<int> const &ix) const
00059   {
00060     size_t addr = 0;
00061     for (size_t i = 0; i < nd; i++) {
00062       addr += ix[i]*static_cast<size_t>(nxc[i]);
00063       if (cvm::debug()) {
00064         if (ix[i] >= nx[i]) {
00065           cvm::error("Error: exceeding bounds in colvar_grid.\n", COLVARS_BUG_ERROR);
00066           return 0;
00067         }
00068       }
00069     }
00070     return addr;
00071   }
00072 
00073 public:
00074 
00076   std::vector<colvarvalue> lower_boundaries;
00077 
00079   std::vector<colvarvalue> upper_boundaries;
00080 
00082   std::vector<bool>        periodic;
00083 
00085   std::vector<bool>        hard_lower_boundaries;
00086 
00088   std::vector<bool>        hard_upper_boundaries;
00089 
00091   std::vector<cvm::real>   widths;
00092 
00094   bool has_parent_data;
00095 
00097   bool has_data;
00098 
00100   inline size_t num_variables() const
00101   {
00102     return nd;
00103   }
00104 
00106   inline std::vector<int> const &number_of_points_vec() const
00107   {
00108     return nx;
00109   }
00110 
00113   inline size_t number_of_points(int const icv = -1) const
00114   {
00115     if (icv < 0) {
00116       return nt;
00117     } else {
00118       return nx[icv];
00119     }
00120   }
00121 
00123   inline std::vector<int> const & sizes() const
00124   {
00125     return nx;
00126   }
00127 
00129   inline void set_sizes(std::vector<int> const &new_sizes)
00130   {
00131     nx = new_sizes;
00132   }
00133 
00135   inline size_t multiplicity() const
00136   {
00137     return mult;
00138   }
00139 
00141   inline void request_actual_value(bool b = true)
00142   {
00143     size_t i;
00144     for (i = 0; i < use_actual_value.size(); i++) {
00145       use_actual_value[i] = b;
00146     }
00147   }
00148 
00150   int setup(std::vector<int> const &nx_i,
00151             T const &t = T(),
00152             size_t const &mult_i = 1)
00153   {
00154     if (cvm::debug()) {
00155       cvm::log("Allocating grid: multiplicity = "+cvm::to_str(mult_i)+
00156                ", dimensions = "+cvm::to_str(nx_i)+".\n");
00157     }
00158 
00159     mult = mult_i;
00160 
00161     data.clear();
00162 
00163     nx = nx_i;
00164     nd = nx.size();
00165 
00166     nxc.resize(nd);
00167 
00168     // setup dimensions
00169     nt = mult;
00170     for (int i = nd-1; i >= 0; i--) {
00171       if (nx[i] <= 0) {
00172         cvm::error("Error: providing an invalid number of grid points, "+
00173                    cvm::to_str(nx[i])+".\n", COLVARS_BUG_ERROR);
00174         return COLVARS_ERROR;
00175       }
00176       nxc[i] = nt;
00177       nt *= nx[i];
00178     }
00179 
00180     if (cvm::debug()) {
00181       cvm::log("Total number of grid elements = "+cvm::to_str(nt)+".\n");
00182     }
00183 
00184     data.reserve(nt);
00185     data.assign(nt, t);
00186 
00187     return COLVARS_OK;
00188   }
00189 
00191   int setup()
00192   {
00193     return setup(this->nx, T(), this->mult);
00194   }
00195 
00197   void reset(T const &t = T())
00198   {
00199     data.assign(nt, t);
00200   }
00201 
00202 
00204   colvar_grid() : has_data(false)
00205   {
00206     nd = nt = 0;
00207     mult = 1;
00208     has_parent_data = false;
00209     this->setup();
00210   }
00211 
00213   virtual ~colvar_grid()
00214   {}
00215 
00219   colvar_grid(colvar_grid<T> const &g) : colvarparse(),
00220                                          nd(g.nd),
00221                                          nx(g.nx),
00222                                          mult(g.mult),
00223                                          data(),
00224                                          cv(g.cv),
00225                                          use_actual_value(g.use_actual_value),
00226                                          lower_boundaries(g.lower_boundaries),
00227                                          upper_boundaries(g.upper_boundaries),
00228                                          periodic(g.periodic),
00229                                          hard_lower_boundaries(g.hard_lower_boundaries),
00230                                          hard_upper_boundaries(g.hard_upper_boundaries),
00231                                          widths(g.widths),
00232                                          has_parent_data(false),
00233                                          has_data(false)
00234   {}
00235 
00240   colvar_grid(std::vector<int> const &nx_i,
00241               T const &t = T(),
00242               size_t mult_i = 1)
00243     : has_parent_data(false), has_data(false)
00244   {
00245     this->setup(nx_i, t, mult_i);
00246   }
00247 
00251   colvar_grid(std::vector<colvar *> const &colvars,
00252               T const &t = T(),
00253               size_t mult_i = 1,
00254               bool add_extra_bin = false)
00255     : has_parent_data(false), has_data(false)
00256   {
00257     (void) t;
00258     this->init_from_colvars(colvars, mult_i, add_extra_bin);
00259   }
00260 
00261   int init_from_colvars(std::vector<colvar *> const &colvars,
00262                         size_t mult_i = 1,
00263                         bool add_extra_bin = false)
00264   {
00265     if (cvm::debug()) {
00266       cvm::log("Reading grid configuration from collective variables.\n");
00267     }
00268 
00269     cv = colvars;
00270     nd = colvars.size();
00271     mult = mult_i;
00272 
00273     size_t i;
00274 
00275     if (cvm::debug()) {
00276       cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+
00277                " collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n");
00278     }
00279 
00280     for (i =  0; i < cv.size(); i++) {
00281 
00282       if (cv[i]->value().type() != colvarvalue::type_scalar) {
00283         cvm::error("Colvar grids can only be automatically "
00284                    "constructed for scalar variables.  "
00285                    "ABF and histogram can not be used; metadynamics "
00286                    "can be used with useGrids disabled.\n", COLVARS_INPUT_ERROR);
00287         return COLVARS_ERROR;
00288       }
00289 
00290       if (cv[i]->width <= 0.0) {
00291         cvm::error("Tried to initialize a grid on a "
00292                    "variable with negative or zero width.\n", COLVARS_INPUT_ERROR);
00293         return COLVARS_ERROR;
00294       }
00295 
00296       widths.push_back(cv[i]->width);
00297       hard_lower_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_lower_boundary));
00298       hard_upper_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_upper_boundary));
00299       periodic.push_back(cv[i]->periodic_boundaries());
00300 
00301       // By default, get reported colvar value (for extended Lagrangian colvars)
00302       use_actual_value.push_back(false);
00303 
00304       // except if a colvar is specified twice in a row
00305       // then the first instance is the actual value
00306       // For histograms of extended-system coordinates
00307       if (i > 0 && cv[i-1] == cv[i]) {
00308         use_actual_value[i-1] = true;
00309       }
00310 
00311       if (add_extra_bin) {
00312         if (periodic[i]) {
00313           // Shift the grid by half the bin width (values at edges instead of center of bins)
00314           lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
00315           upper_boundaries.push_back(cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
00316         } else {
00317           // Make this grid larger by one bin width
00318           lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
00319           upper_boundaries.push_back(cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
00320         }
00321       } else {
00322         lower_boundaries.push_back(cv[i]->lower_boundary);
00323         upper_boundaries.push_back(cv[i]->upper_boundary);
00324       }
00325     }
00326 
00327     this->init_from_boundaries();
00328     return this->setup();
00329   }
00330 
00331   int init_from_boundaries()
00332   {
00333     if (cvm::debug()) {
00334       cvm::log("Configuring grid dimensions from colvars boundaries.\n");
00335     }
00336 
00337     // these will have to be updated
00338     nx.clear();
00339     nxc.clear();
00340     nt = 0;
00341 
00342     for (size_t i =  0; i < lower_boundaries.size(); i++) {
00343       // Re-compute periodicity using current grid boundaries
00344       periodic[i] = cv[i]->periodic_boundaries(lower_boundaries[i].real_value,
00345                                                upper_boundaries[i].real_value);
00346 
00347       cvm::real nbins = ( upper_boundaries[i].real_value -
00348                           lower_boundaries[i].real_value ) / widths[i];
00349       int nbins_round = (int)(nbins+0.5);
00350 
00351       if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
00352         cvm::log("Warning: grid interval("+
00353                  cvm::to_str(lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
00354                  cvm::to_str(upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+
00355                  ") is not commensurate to its bin width("+
00356                  cvm::to_str(widths[i], cvm::cv_width, cvm::cv_prec)+").\n");
00357         upper_boundaries[i].real_value = lower_boundaries[i].real_value +
00358           (nbins_round * widths[i]);
00359       }
00360 
00361       if (cvm::debug())
00362         cvm::log("Number of points is "+cvm::to_str((int) nbins_round)+
00363                  " for the colvar no. "+cvm::to_str(i+1)+".\n");
00364 
00365       nx.push_back(nbins_round);
00366     }
00367 
00368     return COLVARS_OK;
00369   }
00370 
00373   inline void wrap(std::vector<int> & ix) const
00374   {
00375     for (size_t i = 0; i < nd; i++) {
00376       if (periodic[i]) {
00377         ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result
00378       } else {
00379         if (ix[i] < 0 || ix[i] >= nx[i]) {
00380           cvm::error("Trying to wrap illegal index vector (non-PBC) for a grid point: "
00381                      + cvm::to_str(ix), COLVARS_BUG_ERROR);
00382           return;
00383         }
00384       }
00385     }
00386   }
00387 
00390   inline bool wrap_edge(std::vector<int> & ix) const
00391   {
00392     bool edge = false;
00393     for (size_t i = 0; i < nd; i++) {
00394       if (periodic[i]) {
00395         ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result
00396       } else if (ix[i] == -1 || ix[i] == nx[i]) {
00397         edge = true;
00398       }
00399     }
00400     return edge;
00401   }
00402 
00404   inline int current_bin_scalar(int const i) const
00405   {
00406     return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
00407   }
00408 
00411   inline int current_bin_scalar_bound(int const i) const
00412   {
00413     return value_to_bin_scalar_bound(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
00414   }
00415 
00417   inline int current_bin_scalar(int const i, int const iv) const
00418   {
00419     return value_to_bin_scalar(use_actual_value[i] ?
00420                                cv[i]->actual_value().vector1d_value[iv] :
00421                                cv[i]->value().vector1d_value[iv], i);
00422   }
00423 
00426   inline int value_to_bin_scalar(colvarvalue const &value, const int i) const
00427   {
00428     return (int) cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
00429   }
00430 
00432   inline cvm::real current_bin_scalar_fraction(int const i) const
00433   {
00434     return value_to_bin_scalar_fraction(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
00435   }
00436 
00439   inline cvm::real value_to_bin_scalar_fraction(colvarvalue const &value, const int i) const
00440   {
00441     cvm::real x = (value.real_value - lower_boundaries[i].real_value) / widths[i];
00442     return x - cvm::floor(x);
00443   }
00444 
00447   inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
00448   {
00449     int bin_index = cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
00450     if (bin_index < 0) bin_index=0;
00451     if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1;
00452     return (int) bin_index;
00453   }
00454 
00456   inline int value_to_bin_scalar(colvarvalue const &value,
00457                                  colvarvalue const &new_offset,
00458                                  cvm::real   const &new_width) const
00459   {
00460     return (int) cvm::floor( (value.real_value - new_offset.real_value) / new_width );
00461   }
00462 
00465   inline colvarvalue bin_to_value_scalar(int const &i_bin, int const i) const
00466   {
00467     return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
00468   }
00469 
00471   inline colvarvalue bin_to_value_scalar(int const &i_bin,
00472                                          colvarvalue const &new_offset,
00473                                          cvm::real const &new_width) const
00474   {
00475     return new_offset.real_value + new_width * (0.5 + i_bin);
00476   }
00477 
00479   inline void set_value(std::vector<int> const &ix,
00480                         T const &t,
00481                         size_t const &imult = 0)
00482   {
00483     data[this->address(ix)+imult] = t;
00484     has_data = true;
00485   }
00486 
00488   inline void set_value(size_t i, T const &t)
00489   {
00490     data[i] = t;
00491   }
00492 
00497   void delta_grid(colvar_grid<T> const &other_grid)
00498   {
00499 
00500     if (other_grid.multiplicity() != this->multiplicity()) {
00501       cvm::error("Error: trying to subtract two grids with "
00502                  "different multiplicity.\n");
00503       return;
00504     }
00505 
00506     if (other_grid.data.size() != this->data.size()) {
00507       cvm::error("Error: trying to subtract two grids with "
00508                  "different size.\n");
00509       return;
00510     }
00511 
00512     for (size_t i = 0; i < data.size(); i++) {
00513       data[i] = other_grid.data[i] - data[i];
00514     }
00515     has_data = true;
00516   }
00517 
00521   void copy_grid(colvar_grid<T> const &other_grid)
00522   {
00523     if (other_grid.multiplicity() != this->multiplicity()) {
00524       cvm::error("Error: trying to copy two grids with "
00525                  "different multiplicity.\n");
00526       return;
00527     }
00528 
00529     if (other_grid.data.size() != this->data.size()) {
00530       cvm::error("Error: trying to copy two grids with "
00531                  "different size.\n");
00532       return;
00533     }
00534 
00535 
00536     for (size_t i = 0; i < data.size(); i++) {
00537       data[i] = other_grid.data[i];
00538     }
00539     has_data = true;
00540   }
00541 
00544   void raw_data_out(T* out_data) const
00545   {
00546     for (size_t i = 0; i < data.size(); i++) out_data[i] = data[i];
00547   }
00548   void raw_data_out(std::vector<T>& out_data) const
00549   {
00550     out_data = data;
00551   }
00553   void raw_data_in(const T* in_data)
00554   {
00555     for (size_t i = 0; i < data.size(); i++) data[i] = in_data[i];
00556     has_data = true;
00557   }
00558   void raw_data_in(const std::vector<T>& in_data)
00559   {
00560     data = in_data;
00561     has_data = true;
00562   }
00564   size_t raw_data_num() const { return data.size(); }
00565 
00566 
00569   inline T const & value(std::vector<int> const &ix,
00570                          size_t const &imult = 0) const
00571   {
00572     return data[this->address(ix) + imult];
00573   }
00574 
00576   inline T const & value(size_t i) const
00577   {
00578     return data[i];
00579   }
00580 
00582   inline void add_constant(T const &t)
00583   {
00584     for (size_t i = 0; i < nt; i++)
00585       data[i] += t;
00586     has_data = true;
00587   }
00588 
00590   inline void multiply_constant(cvm::real const &a)
00591   {
00592     for (size_t i = 0; i < nt; i++)
00593       data[i] *= a;
00594   }
00595 
00597   inline void remove_small_values(cvm::real const &a)
00598   {
00599     for (size_t i = 0; i < nt; i++)
00600       if(data[i]<a) data[i] = a;
00601   }
00602 
00603 
00606   inline std::vector<int> const get_colvars_index(std::vector<colvarvalue> const &values) const
00607   {
00608     std::vector<int> index = new_index();
00609     for (size_t i = 0; i < nd; i++) {
00610       index[i] = value_to_bin_scalar(values[i], i);
00611     }
00612     return index;
00613   }
00614 
00617   inline std::vector<int> const get_colvars_index() const
00618   {
00619     std::vector<int> index = new_index();
00620     for (size_t i = 0; i < nd; i++) {
00621       index[i] = current_bin_scalar(i);
00622     }
00623     return index;
00624   }
00625 
00628   inline std::vector<int> const get_colvars_index_bound() const
00629   {
00630     std::vector<int> index = new_index();
00631     for (size_t i = 0; i < nd; i++) {
00632       index[i] = current_bin_scalar_bound(i);
00633     }
00634     return index;
00635   }
00636 
00640   inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values,
00641                                                 bool skip_hard_boundaries = false)
00642   {
00643     cvm::real minimum = 1.0E+16;
00644     for (size_t i = 0; i < nd; i++) {
00645 
00646       if (periodic[i]) continue;
00647 
00648       cvm::real dl = cvm::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i];
00649       cvm::real du = cvm::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i];
00650 
00651       if (values[i].real_value < lower_boundaries[i])
00652         dl *= -1.0;
00653       if (values[i].real_value > upper_boundaries[i])
00654         du *= -1.0;
00655 
00656       if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum))
00657         minimum = dl;
00658       if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum))
00659         minimum = du;
00660     }
00661 
00662     return minimum;
00663   }
00664 
00665 
00670   void map_grid(colvar_grid<T> const &other_grid)
00671   {
00672     if (other_grid.multiplicity() != this->multiplicity()) {
00673       cvm::error("Error: trying to merge two grids with values of "
00674                  "different multiplicity.\n");
00675       return;
00676     }
00677 
00678     std::vector<colvarvalue> const &gb  = this->lower_boundaries;
00679     std::vector<cvm::real> const &gw    = this->widths;
00680     std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
00681     std::vector<cvm::real> const &ogw   = other_grid.widths;
00682 
00683     std::vector<int> ix = this->new_index();
00684     std::vector<int> oix = other_grid.new_index();
00685 
00686     if (cvm::debug())
00687       cvm::log("Remapping grid...\n");
00688     for ( ; this->index_ok(ix); this->incr(ix)) {
00689 
00690       for (size_t i = 0; i < nd; i++) {
00691         oix[i] =
00692           value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]),
00693                               ogb[i],
00694                               ogw[i]);
00695       }
00696 
00697       if (! other_grid.index_ok(oix)) {
00698         continue;
00699       }
00700 
00701       for (size_t im = 0; im < mult; im++) {
00702         this->set_value(ix, other_grid.value(oix, im), im);
00703       }
00704     }
00705 
00706     has_data = true;
00707     if (cvm::debug())
00708       cvm::log("Remapping done.\n");
00709   }
00710 
00713   void add_grid(colvar_grid<T> const &other_grid,
00714                 cvm::real scale_factor = 1.0)
00715   {
00716     if (other_grid.multiplicity() != this->multiplicity()) {
00717       cvm::error("Error: trying to sum togetehr two grids with values of "
00718                  "different multiplicity.\n");
00719       return;
00720     }
00721     if (scale_factor != 1.0)
00722       for (size_t i = 0; i < data.size(); i++) {
00723         data[i] += static_cast<T>(scale_factor * other_grid.data[i]);
00724       }
00725     else
00726       // skip multiplication if possible
00727       for (size_t i = 0; i < data.size(); i++) {
00728         data[i] += other_grid.data[i];
00729       }
00730     has_data = true;
00731   }
00732 
00735   virtual T value_output(std::vector<int> const &ix,
00736                          size_t const &imult = 0) const
00737   {
00738     return value(ix, imult);
00739   }
00740 
00744   virtual void value_input(std::vector<int> const &ix,
00745                            T const &t,
00746                            size_t const &imult = 0,
00747                            bool add = false)
00748   {
00749     if ( add )
00750       data[address(ix) + imult] += t;
00751     else
00752       data[address(ix) + imult] = t;
00753     has_data = true;
00754   }
00755 
00756   //   /// Get the pointer to the binned value indexed by ix
00757   //   inline T const *value_p (std::vector<int> const &ix)
00758   //   {
00759   //     return &(data[address (ix)]);
00760   //   }
00761 
00764   inline std::vector<int> const new_index() const
00765   {
00766     return std::vector<int> (nd, 0);
00767   }
00768 
00771   inline bool index_ok(std::vector<int> const &ix) const
00772   {
00773     for (size_t i = 0; i < nd; i++) {
00774       if ( (ix[i] < 0) || (ix[i] >= int(nx[i])) )
00775         return false;
00776     }
00777     return true;
00778   }
00779 
00782   inline void incr(std::vector<int> &ix) const
00783   {
00784     for (int i = ix.size()-1; i >= 0; i--) {
00785 
00786       ix[i]++;
00787 
00788       if (ix[i] >= nx[i]) {
00789 
00790         if (i > 0) {
00791           ix[i] = 0;
00792           continue;
00793         } else {
00794           // this is the last iteration, a non-valid index is being
00795           // set for the outer index, which will be caught by
00796           // index_ok()
00797           ix[0] = nx[0];
00798           return;
00799         }
00800       } else {
00801         return;
00802       }
00803     }
00804   }
00805 
00807   std::ostream & write_params(std::ostream &os)
00808   {
00809     size_t i;
00810     os << "grid_parameters {\n  n_colvars " << nd << "\n";
00811 
00812     os << "  lower_boundaries ";
00813     for (i = 0; i < nd; i++)
00814       os << " " << lower_boundaries[i];
00815     os << "\n";
00816 
00817     os << "  upper_boundaries ";
00818     for (i = 0; i < nd; i++)
00819       os << " " << upper_boundaries[i];
00820     os << "\n";
00821 
00822     os << "  widths ";
00823     for (i = 0; i < nd; i++)
00824       os << " " << widths[i];
00825     os << "\n";
00826 
00827     os << "  sizes ";
00828     for (i = 0; i < nd; i++)
00829       os << " " << nx[i];
00830     os << "\n";
00831 
00832     os << "}\n";
00833     return os;
00834   }
00835 
00837   int parse_params(std::string const &conf,
00838                    colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal)
00839   {
00840     if (cvm::debug()) cvm::log("Reading grid configuration from string.\n");
00841 
00842     std::vector<int> old_nx = nx;
00843     std::vector<colvarvalue> old_lb = lower_boundaries;
00844     std::vector<colvarvalue> old_ub = upper_boundaries;
00845     std::vector<cvm::real> old_w = widths;
00846 
00847     {
00848       size_t nd_in = 0;
00849       // this is only used in state files
00850       colvarparse::get_keyval(conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
00851       if (nd_in != nd) {
00852         cvm::error("Error: trying to read data for a grid "
00853                    "that contains a different number of colvars ("+
00854                    cvm::to_str(nd_in)+") than the grid defined "
00855                    "in the configuration file("+cvm::to_str(nd)+
00856                    ").\n");
00857         return COLVARS_ERROR;
00858       }
00859     }
00860 
00861     // underscore keywords are used in state file
00862     colvarparse::get_keyval(conf, "lower_boundaries",
00863                             lower_boundaries, lower_boundaries, colvarparse::parse_silent);
00864     colvarparse::get_keyval(conf, "upper_boundaries",
00865                             upper_boundaries, upper_boundaries, colvarparse::parse_silent);
00866 
00867     // camel case keywords are used in config file
00868     colvarparse::get_keyval(conf, "lowerBoundaries",
00869                             lower_boundaries, lower_boundaries, parse_mode);
00870     colvarparse::get_keyval(conf, "upperBoundaries",
00871                             upper_boundaries, upper_boundaries, parse_mode);
00872 
00873     colvarparse::get_keyval(conf, "widths", widths, widths, parse_mode);
00874 
00875     // only used in state file
00876     colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent);
00877 
00878     if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
00879 
00880     if (! use_actual_value.size()) use_actual_value.assign(nd, false);
00881     if (! periodic.size()) periodic.assign(nd, false);
00882     if (! widths.size()) widths.assign(nd, 1.0);
00883 
00884     cvm::real eps = 1.e-10;
00885 
00886     bool new_params = false;
00887     if (old_nx.size()) {
00888       for (size_t i = 0; i < nd; i++) {
00889         if (old_nx[i] != nx[i] ||
00890             cvm::sqrt(cv[i]->dist2(old_lb[i], lower_boundaries[i])) > eps ||
00891             cvm::sqrt(cv[i]->dist2(old_ub[i], upper_boundaries[i])) > eps ||
00892             cvm::fabs(old_w[i] - widths[i]) > eps) {
00893           new_params = true;
00894         }
00895       }
00896     } else {
00897       new_params = true;
00898     }
00899 
00900     // reallocate the array in case the grid params have just changed
00901     if (new_params) {
00902       init_from_boundaries();
00903       // data.clear(); // no longer needed: setup calls clear()
00904       return this->setup(nx, T(), mult);
00905     }
00906 
00907     return COLVARS_OK;
00908   }
00909 
00913   void check_consistency()
00914   {
00915     for (size_t i = 0; i < nd; i++) {
00916       if ( (cvm::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
00917                                    lower_boundaries[i])) > 1.0E-10) ||
00918            (cvm::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
00919                                    upper_boundaries[i])) > 1.0E-10) ||
00920            (cvm::sqrt(cv[i]->dist2(cv[i]->width,
00921                                    widths[i])) > 1.0E-10) ) {
00922         cvm::error("Error: restart information for a grid is "
00923                    "inconsistent with that of its colvars.\n");
00924         return;
00925       }
00926     }
00927   }
00928 
00929 
00932   void check_consistency(colvar_grid<T> const &other_grid)
00933   {
00934     for (size_t i = 0; i < nd; i++) {
00935       // we skip dist2(), because periodicities and the like should
00936       // matter: boundaries should be EXACTLY the same (otherwise,
00937       // map_grid() should be used)
00938       if ( (cvm::fabs(other_grid.lower_boundaries[i] -
00939                       lower_boundaries[i]) > 1.0E-10) ||
00940            (cvm::fabs(other_grid.upper_boundaries[i] -
00941                       upper_boundaries[i]) > 1.0E-10) ||
00942            (cvm::fabs(other_grid.widths[i] -
00943                       widths[i]) > 1.0E-10) ||
00944            (data.size() != other_grid.data.size()) ) {
00945         cvm::error("Error: inconsistency between "
00946                    "two grids that are supposed to be equal, "
00947                    "aside from the data stored.\n");
00948         return;
00949       }
00950     }
00951   }
00952 
00953 
00955   std::istream & read_restart(std::istream &is)
00956   {
00957     std::streampos const start_pos = is.tellg();
00958     std::string key, conf;
00959     if ((is >> key) && (key == std::string("grid_parameters"))) {
00960       is.seekg(start_pos, std::ios::beg);
00961       is >> colvarparse::read_block("grid_parameters", &conf);
00962       parse_params(conf, colvarparse::parse_silent);
00963     } else {
00964       cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
00965       is.seekg(start_pos, std::ios::beg);
00966     }
00967     read_raw(is);
00968     return is;
00969   }
00970 
00972   std::ostream & write_restart(std::ostream &os)
00973   {
00974     write_params(os);
00975     write_raw(os);
00976     return os;
00977   }
00978 
00979 
00983   std::ostream & write_raw(std::ostream &os,
00984                            size_t const buf_size = 3) const
00985   {
00986     std::streamsize const w = os.width();
00987     std::streamsize const p = os.precision();
00988 
00989     std::vector<int> ix = new_index();
00990     size_t count = 0;
00991     for ( ; index_ok(ix); incr(ix)) {
00992       for (size_t imult = 0; imult < mult; imult++) {
00993         os << " "
00994            << std::setw(w) << std::setprecision(p)
00995            << value_output(ix, imult);
00996         if (((++count) % buf_size) == 0)
00997           os << "\n";
00998       }
00999     }
01000     // write a final newline only if buffer is not empty
01001     if ((count % buf_size) != 0)
01002       os << "\n";
01003 
01004     return os;
01005   }
01006 
01008   std::istream & read_raw(std::istream &is)
01009   {
01010     std::streampos const start_pos = is.tellg();
01011 
01012     for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
01013       for (size_t imult = 0; imult < mult; imult++) {
01014         T new_value;
01015         if (is >> new_value) {
01016           value_input(ix, new_value, imult);
01017         } else {
01018           is.clear();
01019           is.seekg(start_pos, std::ios::beg);
01020           is.setstate(std::ios::failbit);
01021           cvm::error("Error: failed to read all of the grid points from file.  Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n");
01022           return is;
01023         }
01024       }
01025     }
01026 
01027     has_data = true;
01028     return is;
01029   }
01030 
01032   std::istream & read_multicol(std::istream &is, bool add = false);
01033 
01035   int read_multicol(std::string const &filename,
01036                     std::string description = "grid file",
01037                     bool add = false);
01038 
01040   std::ostream & write_multicol(std::ostream &os) const;
01041 
01043   int write_multicol(std::string const &filename,
01044                      std::string description = "grid file") const;
01045 
01047   std::ostream & write_opendx(std::ostream &os) const;
01048 
01050   int write_opendx(std::string const &filename,
01051                    std::string description = "grid file") const;
01052 };
01053 
01054 
01055 
01058 class colvar_grid_count : public colvar_grid<size_t>
01059 {
01060 public:
01061 
01063   colvar_grid_count();
01064 
01066   virtual ~colvar_grid_count()
01067   {}
01068 
01070   colvar_grid_count(std::vector<int> const &nx_i,
01071                     size_t const           &def_count = 0);
01072 
01074   colvar_grid_count(std::vector<colvar *>  &colvars,
01075                     size_t const           &def_count = 0,
01076                     bool                   add_extra_bin = false);
01077 
01079   inline void incr_count(std::vector<int> const &ix)
01080   {
01081     ++(data[this->address(ix)]);
01082   }
01083 
01085   inline size_t const & new_count(std::vector<int> const &ix,
01086                                   size_t const &imult = 0)
01087   {
01088     return new_data[address(ix) + imult];
01089   }
01090 
01092   std::istream & read_multicol(std::istream &is, bool add = false);
01093 
01095   int read_multicol(std::string const &filename,
01096                             std::string description = "grid file",
01097                             bool add = false);
01098 
01100   std::ostream & write_multicol(std::ostream &os) const;
01101 
01103   int write_multicol(std::string const &filename,
01104                      std::string description = "grid file") const;
01105 
01107   std::ostream & write_opendx(std::ostream &os) const;
01108 
01110   int write_opendx(std::string const &filename,
01111                    std::string description = "grid file") const;
01112 
01114   virtual void value_input(std::vector<int> const &ix,
01115                            size_t const &t,
01116                            size_t const &imult = 0,
01117                            bool add = false)
01118   {
01119     (void) imult;
01120     if (add) {
01121       data[address(ix)] += t;
01122       if (this->has_parent_data) {
01123         // save newly read data for inputting parent grid
01124         new_data[address(ix)] = t;
01125       }
01126     } else {
01127       data[address(ix)] = t;
01128     }
01129     has_data = true;
01130   }
01131 
01134   inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
01135                                             int n = 0)
01136   {
01137     cvm::real A0, A1, A2;
01138     std::vector<int> ix = ix0;
01139 
01140     // TODO this can be rewritten more concisely with wrap_edge()
01141     if (periodic[n]) {
01142       ix[n]--; wrap(ix);
01143       A0 = value(ix);
01144       ix = ix0;
01145       ix[n]++; wrap(ix);
01146       A1 = value(ix);
01147       if (A0 * A1 == 0) {
01148         return 0.; // can't handle empty bins
01149       } else {
01150         return (cvm::logn(A1) - cvm::logn(A0))
01151           / (widths[n] * 2.);
01152       }
01153     } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
01154       ix[n]--;
01155       A0 = value(ix);
01156       ix = ix0;
01157       ix[n]++;
01158       A1 = value(ix);
01159       if (A0 * A1 == 0) {
01160         return 0.; // can't handle empty bins
01161       } else {
01162         return (cvm::logn(A1) - cvm::logn(A0))
01163           / (widths[n] * 2.);
01164       }
01165     } else {
01166       // edge: use 2nd order derivative
01167       int increment = (ix[n] == 0 ? 1 : -1);
01168       // move right from left edge, or the other way around
01169       A0 = value(ix);
01170       ix[n] += increment; A1 = value(ix);
01171       ix[n] += increment; A2 = value(ix);
01172       if (A0 * A1 * A2 == 0) {
01173         return 0.; // can't handle empty bins
01174       } else {
01175         return (-1.5 * cvm::logn(A0) + 2. * cvm::logn(A1)
01176           - 0.5 * cvm::logn(A2)) * increment / widths[n];
01177       }
01178     }
01179   }
01180 
01183   inline cvm::real gradient_finite_diff(const std::vector<int> &ix0,
01184                                         int n = 0)
01185   {
01186     cvm::real A0, A1, A2;
01187     std::vector<int> ix = ix0;
01188 
01189     // FIXME this can be rewritten more concisely with wrap_edge()
01190     if (periodic[n]) {
01191       ix[n]--; wrap(ix);
01192       A0 = value(ix);
01193       ix = ix0;
01194       ix[n]++; wrap(ix);
01195       A1 = value(ix);
01196       if (A0 * A1 == 0) {
01197         return 0.; // can't handle empty bins
01198       } else {
01199         return (A1 - A0) / (widths[n] * 2.);
01200       }
01201     } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
01202       ix[n]--;
01203       A0 = value(ix);
01204       ix = ix0;
01205       ix[n]++;
01206       A1 = value(ix);
01207       if (A0 * A1 == 0) {
01208         return 0.; // can't handle empty bins
01209       } else {
01210         return (A1 - A0) / (widths[n] * 2.);
01211       }
01212     } else {
01213       // edge: use 2nd order derivative
01214       int increment = (ix[n] == 0 ? 1 : -1);
01215       // move right from left edge, or the other way around
01216       A0 = value(ix);
01217       ix[n] += increment; A1 = value(ix);
01218       ix[n] += increment; A2 = value(ix);
01219       return (-1.5 * A0 + 2. * A1
01220           - 0.5 * A2) * increment / widths[n];
01221     }
01222   }
01223 };
01224 
01225 
01227 class colvar_grid_scalar : public colvar_grid<cvm::real>
01228 {
01229 public:
01230 
01233   colvar_grid_count *samples;
01234 
01236   colvar_grid_scalar();
01237 
01239   colvar_grid_scalar(colvar_grid_scalar const &g);
01240 
01242   virtual ~colvar_grid_scalar();
01243 
01245   colvar_grid_scalar(std::vector<int> const &nx_i);
01246 
01248   colvar_grid_scalar(std::vector<colvar *> &colvars,
01249                      bool add_extra_bin = false);
01250 
01252   inline void acc_value(std::vector<int> const &ix,
01253                         cvm::real const &new_value,
01254                         size_t const &imult = 0)
01255   {
01256     (void) imult;
01257     // only legal value of imult here is 0
01258     data[address(ix)] += new_value;
01259     if (samples)
01260       samples->incr_count(ix);
01261     has_data = true;
01262   }
01263 
01265   std::istream & read_multicol(std::istream &is, bool add = false);
01266 
01268   int read_multicol(std::string const &filename,
01269                     std::string description = "grid file",
01270                     bool add = false);
01271 
01273   std::ostream & write_multicol(std::ostream &os) const;
01274 
01276   int write_multicol(std::string const &filename,
01277                      std::string description = "grid file") const;
01278 
01280   std::ostream & write_opendx(std::ostream &os) const;
01281 
01283   int write_opendx(std::string const &filename,
01284                    std::string description = "grid file") const;
01285 
01290   inline void vector_gradient_finite_diff( const std::vector<int> &ix0, std::vector<cvm::real> &grad)
01291   {
01292     cvm::real A0, A1;
01293     std::vector<int> ix;
01294     size_t i, j, k, n;
01295 
01296     if (nd == 2) {
01297       for (n = 0; n < 2; n++) {
01298         ix = ix0;
01299         A0 = value(ix);
01300         ix[n]++; wrap(ix);
01301         A1 = value(ix);
01302         ix[1-n]++; wrap(ix);
01303         A1 += value(ix);
01304         ix[n]--; wrap(ix);
01305         A0 += value(ix);
01306         grad[n] = 0.5 * (A1 - A0) / widths[n];
01307       }
01308     } else if (nd == 3) {
01309 
01310       cvm::real p[8]; // potential values within cube, indexed in binary (4 i + 2 j + k)
01311       ix = ix0;
01312       int index = 0;
01313       for (i = 0; i<2; i++) {
01314         ix[1] = ix0[1];
01315         for (j = 0; j<2; j++) {
01316           ix[2] = ix0[2];
01317           for (k = 0; k<2; k++) {
01318             wrap(ix);
01319             p[index++] = value(ix);
01320             ix[2]++;
01321           }
01322           ix[1]++;
01323         }
01324         ix[0]++;
01325       }
01326 
01327       // The following would be easier to read using binary literals
01328       //                  100    101    110    111      000    001    010   011
01329       grad[0] = 0.25 * ((p[4] + p[5] + p[6] + p[7]) - (p[0] + p[1] + p[2] + p[3])) / widths[0];
01330       //                  010     011    110   111      000    001    100   101
01331       grad[1] = 0.25 * ((p[2] + p[3] + p[6] + p[7]) - (p[0] + p[1] + p[4] + p[5])) / widths[1];
01332       //                  001    011     101   111      000    010   100    110
01333       grad[2] = 0.25 * ((p[1] + p[3] + p[5] + p[7]) - (p[0] + p[2] + p[4] + p[6])) / widths[2];
01334     } else {
01335       cvm::error("Finite differences available in dimension 2 and 3 only.");
01336     }
01337   }
01338 
01341   inline cvm::real gradient_finite_diff(const std::vector<int> &ix0,
01342                                         int n = 0)
01343   {
01344     cvm::real A0, A1, A2;
01345     std::vector<int> ix = ix0;
01346 
01347     // FIXME this can be rewritten more concisely with wrap_edge()
01348     if (periodic[n]) {
01349       ix[n]--; wrap(ix);
01350       A0 = value(ix);
01351       ix = ix0;
01352       ix[n]++; wrap(ix);
01353       A1 = value(ix);
01354       if (A0 * A1 == 0) {
01355         return 0.; // can't handle empty bins
01356       } else {
01357         return (A1 - A0) / (widths[n] * 2.);
01358       }
01359     } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
01360       ix[n]--;
01361       A0 = value(ix);
01362       ix = ix0;
01363       ix[n]++;
01364       A1 = value(ix);
01365       if (A0 * A1 == 0) {
01366         return 0.; // can't handle empty bins
01367       } else {
01368         return (A1 - A0) / (widths[n] * 2.);
01369       }
01370     } else {
01371       // edge: use 2nd order derivative
01372       int increment = (ix[n] == 0 ? 1 : -1);
01373       // move right from left edge, or the other way around
01374       A0 = value(ix);
01375       ix[n] += increment; A1 = value(ix);
01376       ix[n] += increment; A2 = value(ix);
01377       return (-1.5 * A0 + 2. * A1
01378           - 0.5 * A2) * increment / widths[n];
01379     }
01380   }
01381 
01384   virtual cvm::real value_output(std::vector<int> const &ix,
01385                                  size_t const &imult = 0) const
01386   {
01387     if (imult > 0) {
01388       cvm::error("Error: trying to access a component "
01389                  "larger than 1 in a scalar data grid.\n");
01390       return 0.;
01391     }
01392     if (samples) {
01393       return (samples->value(ix) > 0) ?
01394         (data[address(ix)] / cvm::real(samples->value(ix))) :
01395         0.0;
01396     } else {
01397       return data[address(ix)];
01398     }
01399   }
01400 
01402   virtual void value_input(std::vector<int> const &ix,
01403                            cvm::real const &new_value,
01404                            size_t const &imult = 0,
01405                            bool add = false)
01406   {
01407     if (imult > 0) {
01408       cvm::error("Error: trying to access a component "
01409                  "larger than 1 in a scalar data grid.\n");
01410       return;
01411     }
01412     if (add) {
01413       if (samples)
01414         data[address(ix)] += new_value * samples->new_count(ix);
01415       else
01416         data[address(ix)] += new_value;
01417     } else {
01418       if (samples)
01419         data[address(ix)] = new_value * samples->value(ix);
01420       else
01421         data[address(ix)] = new_value;
01422     }
01423     has_data = true;
01424   }
01425 
01427   cvm::real maximum_value() const;
01428 
01430   cvm::real minimum_value() const;
01431 
01433   cvm::real minimum_pos_value() const;
01434 
01436   cvm::real integral() const;
01437 
01440   cvm::real entropy() const;
01441 };
01442 
01443 
01444 
01446 class colvar_grid_gradient : public colvar_grid<cvm::real>
01447 {
01448 public:
01449 
01452   colvar_grid_count *samples;
01453 
01456   colvar_grid_scalar *weights;
01457 
01459   colvar_grid_gradient();
01460 
01462   virtual ~colvar_grid_gradient()
01463   {}
01464 
01466   colvar_grid_gradient(std::vector<int> const &nx_i);
01467 
01469   colvar_grid_gradient(std::vector<colvar *>  &colvars);
01470 
01472   colvar_grid_gradient(std::string &filename);
01473 
01475   virtual std::istream & read_multicol(std::istream &is, bool add = false);
01476 
01478   virtual int read_multicol(std::string const &filename,
01479                             std::string description = "grid file",
01480                             bool add = false);
01481 
01483   virtual std::ostream & write_multicol(std::ostream &os) const;
01484 
01486   virtual int write_multicol(std::string const &filename,
01487                              std::string description = "grid file") const;
01488 
01490   virtual std::ostream & write_opendx(std::ostream &os) const;
01491 
01493   virtual int write_opendx(std::string const &filename,
01494                            std::string description = "grid file") const;
01495 
01497   inline void vector_value(std::vector<int> const &ix, std::vector<cvm::real> &v) const
01498   {
01499     cvm::real const * p = &value(ix);
01500     if (samples) {
01501       int count = samples->value(ix);
01502       if (count) {
01503         cvm::real invcount = 1.0 / count;
01504         for (size_t i = 0; i < mult; i++) {
01505           v[i] = invcount * p[i];
01506         }
01507       } else {
01508         for (size_t i = 0; i < mult; i++) {
01509           v[i] = 0.0;
01510         }
01511       }
01512     } else {
01513       for (size_t i = 0; i < mult; i++) {
01514         v[i] = p[i];
01515       }
01516     }
01517   }
01518 
01520   inline void acc_value(std::vector<int> const &ix, std::vector<colvarvalue> const &values) {
01521     for (size_t imult = 0; imult < mult; imult++) {
01522       data[address(ix) + imult] += values[imult].real_value;
01523     }
01524     if (samples)
01525       samples->incr_count(ix);
01526   }
01527 
01530   inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) {
01531     for (size_t imult = 0; imult < mult; imult++) {
01532       data[address(ix) + imult] -= forces[imult];
01533     }
01534     if (samples)
01535       samples->incr_count(ix);
01536   }
01537 
01540   inline void acc_force_weighted(std::vector<int> const &ix,
01541                                  cvm::real const *forces,
01542                                  cvm::real weight) {
01543     for (size_t imult = 0; imult < mult; imult++) {
01544       data[address(ix) + imult] -= forces[imult] * weight;
01545     }
01546     weights->acc_value(ix, weight);
01547   }
01548 
01551   virtual cvm::real value_output(std::vector<int> const &ix,
01552                                  size_t const &imult = 0) const
01553   {
01554     if (samples)
01555       return (samples->value(ix) > 0) ?
01556         (data[address(ix) + imult] / cvm::real(samples->value(ix))) :
01557         0.0;
01558     else
01559       return data[address(ix) + imult];
01560   }
01561 
01565   virtual void value_input(std::vector<int> const &ix,
01566                            cvm::real const &new_value,
01567                            size_t const &imult = 0,
01568                            bool add = false)
01569   {
01570     if (add) {
01571       if (samples)
01572         data[address(ix) + imult] += new_value * samples->new_count(ix);
01573       else
01574         data[address(ix) + imult] += new_value;
01575     } else {
01576       if (samples)
01577         data[address(ix) + imult] = new_value * samples->value(ix);
01578       else
01579         data[address(ix) + imult] = new_value;
01580     }
01581     has_data = true;
01582   }
01583 
01584 
01586   inline cvm::real average()
01587   {
01588     size_t n = 0;
01589 
01590     if (nd != 1 || nx[0] == 0) {
01591       return 0.0;
01592     }
01593 
01594     cvm::real sum = 0.0;
01595     std::vector<int> ix = new_index();
01596     if (samples) {
01597       for ( ; index_ok(ix); incr(ix)) {
01598         if ( (n = samples->value(ix)) )
01599           sum += value(ix) / n;
01600       }
01601     } else {
01602       for ( ; index_ok(ix); incr(ix)) {
01603         sum += value(ix);
01604       }
01605     }
01606     return (sum / cvm::real(nx[0]));
01607   }
01608 
01611   void write_1D_integral(std::ostream &os);
01612 
01613 };
01614 
01615 
01616 
01618 
01619 class integrate_potential : public colvar_grid_scalar
01620 {
01621   public:
01622 
01623   integrate_potential();
01624 
01625   virtual ~integrate_potential()
01626   {}
01627 
01629   integrate_potential(std::vector<colvar *> &colvars, colvar_grid_gradient * gradients);
01630 
01632   integrate_potential(colvar_grid_gradient * gradients);
01633 
01635   int integrate(const int itmax, const cvm::real & tol, cvm::real & err);
01636 
01639   void update_div_neighbors(const std::vector<int> &ix);
01640 
01643   void set_div();
01644 
01647   inline void set_zero_minimum() {
01648     add_constant(-1.0 * minimum_value());
01649   }
01650 
01651   protected:
01652 
01653   // Reference to gradient grid
01654   colvar_grid_gradient *gradients;
01655 
01657   std::vector<cvm::real> divergence;
01658 
01659 //   std::vector<cvm::real> inv_lap_diag; // Inverse of the diagonal of the Laplacian; for conditioning
01660 
01663   void update_div_local(const std::vector<int> &ix);
01664 
01668   void get_grad(cvm::real * g, std::vector<int> &ix);
01669 
01671   void nr_linbcg_sym(const std::vector<cvm::real> &b, std::vector<cvm::real> &x,
01672                      const cvm::real &tol, const int itmax, int &iter, cvm::real &err);
01673 
01675   cvm::real l2norm(const std::vector<cvm::real> &x);
01676 
01678   void atimes(const std::vector<cvm::real> &x, std::vector<cvm::real> &r);
01679 
01680 //   /// Inversion of preconditioner matrix
01681 //   void asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x);
01682 };
01683 
01684 #endif
01685 

Generated on Thu Mar 28 02:42:54 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002