colvargrid.C

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 <ctime>
00011 #include <fstream>
00012 
00013 #include "colvarmodule.h"
00014 #include "colvarvalue.h"
00015 #include "colvarparse.h"
00016 #include "colvar.h"
00017 #include "colvarcomp.h"
00018 #include "colvargrid.h"
00019 #include "colvargrid_def.h"
00020 
00021 
00022 
00023 colvar_grid_count::colvar_grid_count()
00024 : colvar_grid<size_t>()
00025 {
00026 mult = 1;
00027 }
00028 
00029 colvar_grid_count::colvar_grid_count(std::vector<int> const &nx_i,
00030 size_t const &def_count)
00031 : colvar_grid<size_t>(nx_i, def_count, 1)
00032 {}
00033 
00034 colvar_grid_count::colvar_grid_count(std::vector<colvar *> &colvars,
00035 size_t const &def_count,
00036 bool margin)
00037 : colvar_grid<size_t>(colvars, def_count, 1, margin)
00038 {}
00039 
00040 std::istream & colvar_grid_count::read_multicol(std::istream &is, bool add)
00041 {
00042 return colvar_grid<size_t>::read_multicol(is, add);
00043 }
00044 
00045 int colvar_grid_count::read_multicol(std::string const &filename,
00046 std::string description,
00047 bool add)
00048 {
00049 return colvar_grid<size_t>::read_multicol(filename, description, add);
00050 }
00051 
00052 std::ostream & colvar_grid_count::write_multicol(std::ostream &os) const
00053 {
00054 return colvar_grid<size_t>::write_multicol(os);
00055 }
00056 
00057 int colvar_grid_count::write_multicol(std::string const &filename,
00058 std::string description) const
00059 {
00060 return colvar_grid<size_t>::write_multicol(filename, description);
00061 }
00062 
00063 std::ostream & colvar_grid_count::write_opendx(std::ostream &os) const
00064 {
00065 return colvar_grid<size_t>::write_opendx(os);
00066 }
00067 
00068 int colvar_grid_count::write_opendx(std::string const &filename,
00069 std::string description) const
00070 {
00071 return colvar_grid<size_t>::write_opendx(filename, description);
00072 }
00073 
00074 
00075 
00076 colvar_grid_scalar::colvar_grid_scalar()
00077 : colvar_grid<cvm::real>(), samples(NULL)
00078 {}
00079 
00080 colvar_grid_scalar::colvar_grid_scalar(colvar_grid_scalar const &g)
00081 : colvar_grid<cvm::real>(g), samples(NULL)
00082 {
00083 }
00084 
00085 colvar_grid_scalar::colvar_grid_scalar(std::vector<int> const &nx_i)
00086 : colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL)
00087 {
00088 }
00089 
00090 colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool margin)
00091 : colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL)
00092 {
00093 }
00094 
00095 colvar_grid_scalar::~colvar_grid_scalar()
00096 {
00097 }
00098 
00099 std::istream & colvar_grid_scalar::read_multicol(std::istream &is, bool add)
00100 {
00101 return colvar_grid<cvm::real>::read_multicol(is, add);
00102 }
00103 
00104 int colvar_grid_scalar::read_multicol(std::string const &filename,
00105 std::string description,
00106 bool add)
00107 {
00108 return colvar_grid<cvm::real>::read_multicol(filename, description, add);
00109 }
00110 
00111 std::ostream & colvar_grid_scalar::write_multicol(std::ostream &os) const
00112 {
00113 return colvar_grid<cvm::real>::write_multicol(os);
00114 }
00115 
00116 int colvar_grid_scalar::write_multicol(std::string const &filename,
00117 std::string description) const
00118 {
00119 return colvar_grid<cvm::real>::write_multicol(filename, description);
00120 }
00121 
00122 std::ostream & colvar_grid_scalar::write_opendx(std::ostream &os) const
00123 {
00124 return colvar_grid<cvm::real>::write_opendx(os);
00125 }
00126 
00127 int colvar_grid_scalar::write_opendx(std::string const &filename,
00128 std::string description) const
00129 {
00130 return colvar_grid<cvm::real>::write_opendx(filename, description);
00131 }
00132 
00133 
00134 cvm::real colvar_grid_scalar::maximum_value() const
00135 {
00136 cvm::real max = data[0];
00137 for (size_t i = 0; i < nt; i++) {
00138 if (data[i] > max) max = data[i];
00139 }
00140 return max;
00141 }
00142 
00143 
00144 cvm::real colvar_grid_scalar::minimum_value() const
00145 {
00146 cvm::real min = data[0];
00147 for (size_t i = 0; i < nt; i++) {
00148 if (data[i] < min) min = data[i];
00149 }
00150 return min;
00151 }
00152 
00153 cvm::real colvar_grid_scalar::minimum_pos_value() const
00154 {
00155 cvm::real minpos = data[0];
00156 size_t i;
00157 for (i = 0; i < nt; i++) {
00158 if(data[i] > 0) {
00159 minpos = data[i];
00160 break;
00161 }
00162 }
00163 for (i = 0; i < nt; i++) {
00164 if (data[i] > 0 && data[i] < minpos) minpos = data[i];
00165 }
00166 return minpos;
00167 }
00168 
00169 cvm::real colvar_grid_scalar::integral() const
00170 {
00171 cvm::real sum = 0.0;
00172 for (size_t i = 0; i < nt; i++) {
00173 sum += data[i];
00174 }
00175 cvm::real bin_volume = 1.0;
00176 for (size_t id = 0; id < widths.size(); id++) {
00177 bin_volume *= widths[id];
00178 }
00179 return bin_volume * sum;
00180 }
00181 
00182 
00183 cvm::real colvar_grid_scalar::entropy() const
00184 {
00185 cvm::real sum = 0.0;
00186 for (size_t i = 0; i < nt; i++) {
00187 if (data[i] >0) {
00188 sum += -1.0 * data[i] * cvm::logn(data[i]);
00189 }
00190 }
00191 cvm::real bin_volume = 1.0;
00192 for (size_t id = 0; id < widths.size(); id++) {
00193 bin_volume *= widths[id];
00194 }
00195 return bin_volume * sum;
00196 }
00197 
00198 
00199 colvar_grid_gradient::colvar_grid_gradient()
00200 : colvar_grid<cvm::real>(),
00201 samples(NULL),
00202 weights(NULL)
00203 {}
00204 
00205 colvar_grid_gradient::colvar_grid_gradient(std::vector<int> const &nx_i)
00206 : colvar_grid<cvm::real>(nx_i, 0.0, nx_i.size()),
00207 samples(NULL),
00208 weights(NULL)
00209 {}
00210 
00211 colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars)
00212 : colvar_grid<cvm::real>(colvars, 0.0, colvars.size()),
00213 samples(NULL),
00214 weights(NULL)
00215 {}
00216 
00217 
00218 colvar_grid_gradient::colvar_grid_gradient(std::string &filename)
00219 : colvar_grid<cvm::real>(),
00220 samples(NULL),
00221 weights(NULL)
00222 {
00223 std::istream &is = cvm::main()->proxy->input_stream(filename,
00224 "gradient file");
00225 if (!is) {
00226 return;
00227 }
00228 
00229 // Data in the header: nColvars, then for each
00230 // xiMin, dXi, nPoints, periodic flag
00231 
00232 std::string hash;
00233 size_t i;
00234 
00235 if ( !(is >> hash) || (hash != "#") ) {
00236 cvm::error("Error reading grid at position "+
00237 cvm::to_str(static_cast<size_t>(is.tellg()))+
00238 " in stream(read \"" + hash + "\")\n");
00239 return;
00240 }
00241 
00242 is >> nd;
00243 
00244 if (nd > 50) {
00245 cvm::error("Error: excessive number of dimensions in file \""+
00246 filename+"\". Please ensure that the file is not corrupt.\n",
00247 COLVARS_INPUT_ERROR);
00248 return;
00249 }
00250 
00251 mult = nd;
00252 std::vector<cvm::real> lower_in(nd), widths_in(nd);
00253 std::vector<int> nx_in(nd);
00254 std::vector<int> periodic_in(nd);
00255 
00256 for (i = 0; i < nd; i++ ) {
00257 if ( !(is >> hash) || (hash != "#") ) {
00258 cvm::error("Error reading grid at position "+
00259 cvm::to_str(static_cast<size_t>(is.tellg()))+
00260 " in stream(read \"" + hash + "\")\n");
00261 return;
00262 }
00263 
00264 is >> lower_in[i] >> widths_in[i] >> nx_in[i] >> periodic_in[i];
00265 }
00266 
00267 this->setup(nx_in, 0., mult);
00268 
00269 widths = widths_in;
00270 
00271 for (i = 0; i < nd; i++ ) {
00272 lower_boundaries.push_back(colvarvalue(lower_in[i]));
00273 periodic.push_back(static_cast<bool>(periodic_in[i]));
00274 }
00275 
00276 // Reset the istream for read_multicol, which expects the whole file
00277 is.clear();
00278 is.seekg(0);
00279 read_multicol(is);
00280 cvm::main()->proxy->close_input_stream(filename);
00281 }
00282 
00283 
00284 std::istream & colvar_grid_gradient::read_multicol(std::istream &is, bool add)
00285 {
00286 return colvar_grid<cvm::real>::read_multicol(is, add);
00287 }
00288 
00289 int colvar_grid_gradient::read_multicol(std::string const &filename,
00290 std::string description,
00291 bool add)
00292 {
00293 return colvar_grid<cvm::real>::read_multicol(filename, description, add);
00294 }
00295 
00296 std::ostream & colvar_grid_gradient::write_multicol(std::ostream &os) const
00297 {
00298 return colvar_grid<cvm::real>::write_multicol(os);
00299 }
00300 
00301 int colvar_grid_gradient::write_multicol(std::string const &filename,
00302 std::string description) const
00303 {
00304 return colvar_grid<cvm::real>::write_multicol(filename, description);
00305 }
00306 
00307 std::ostream & colvar_grid_gradient::write_opendx(std::ostream &os) const
00308 {
00309 return colvar_grid<cvm::real>::write_opendx(os);
00310 }
00311 
00312 int colvar_grid_gradient::write_opendx(std::string const &filename,
00313 std::string description) const
00314 {
00315 return colvar_grid<cvm::real>::write_opendx(filename, description);
00316 }
00317 
00318 
00319 void colvar_grid_gradient::write_1D_integral(std::ostream &os)
00320 {
00321 cvm::real bin, min, integral;
00322 std::vector<cvm::real> int_vals;
00323 
00324 os << "# xi A(xi)\n";
00325 
00326 if (cv.size() != 1) {
00327 cvm::error("Cannot write integral for multi-dimensional gradient grids.");
00328 return;
00329 }
00330 
00331 integral = 0.0;
00332 int_vals.push_back(0.0);
00333 min = 0.0;
00334 
00335 // correction for periodic colvars, so that the PMF is periodic
00336 cvm::real corr;
00337 if (periodic[0]) {
00338 corr = average();
00339 } else {
00340 corr = 0.0;
00341 }
00342 
00343 for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
00344 
00345 if (samples) {
00346 size_t const samples_here = samples->value(ix);
00347 if (samples_here)
00348 integral += (value(ix) / cvm::real(samples_here) - corr) * cv[0]->width;
00349 } else {
00350 integral += (value(ix) - corr) * cv[0]->width;
00351 }
00352 
00353 if ( integral < min ) min = integral;
00354 int_vals.push_back(integral);
00355 }
00356 
00357 bin = 0.0;
00358 for ( int i = 0; i < nx[0]; i++, bin += 1.0 ) {
00359 os << std::setw(10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
00360 << std::setw(cvm::cv_width)
00361 << std::setprecision(cvm::cv_prec)
00362 << int_vals[i] - min << "\n";
00363 }
00364 
00365 os << std::setw(10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
00366 << std::setw(cvm::cv_width)
00367 << std::setprecision(cvm::cv_prec)
00368 << int_vals[nx[0]] - min << "\n";
00369 
00370 return;
00371 }
00372 
00373 
00374 
00375 integrate_potential::integrate_potential(std::vector<colvar *> &colvars, colvar_grid_gradient * gradients)
00376 : colvar_grid_scalar(colvars, true),
00377 gradients(gradients)
00378 {
00379 // parent class colvar_grid_scalar is constructed with margin option set to true
00380 // hence PMF grid is wider than gradient grid if non-PBC
00381 
00382 if (nd > 1) {
00383 cvm::main()->cite_feature("Poisson integration of 2D/3D free energy surfaces");
00384 divergence.resize(nt);
00385 
00386 // Compute inverse of Laplacian diagonal for Jacobi preconditioning
00387 // For now all code related to preconditioning is commented out
00388 // until a method better than Jacobi is implemented
00389 // cvm::log("Preparing inverse diagonal for preconditioning...\n");
00390 // inv_lap_diag.resize(nt);
00391 // std::vector<cvm::real> id(nt), lap_col(nt);
00392 // for (int i = 0; i < nt; i++) {
00393 // if (i % (nt / 100) == 0)
00394 // cvm::log(cvm::to_str(i));
00395 // id[i] = 1.;
00396 // atimes(id, lap_col);
00397 // id[i] = 0.;
00398 // inv_lap_diag[i] = 1. / lap_col[i];
00399 // }
00400 // cvm::log("Done.\n");
00401 }
00402 }
00403 
00404 
00405 integrate_potential::integrate_potential(colvar_grid_gradient * gradients)
00406 : gradients(gradients)
00407 {
00408 nd = gradients->num_variables();
00409 nx = gradients->number_of_points_vec();
00410 widths = gradients->widths;
00411 periodic = gradients->periodic;
00412 
00413 // Expand grid by 1 bin in non-periodic dimensions
00414 for (size_t i = 0; i < nd; i++ ) {
00415 if (!periodic[i]) nx[i]++;
00416 // Shift the grid by half the bin width (values at edges instead of center of bins)
00417 lower_boundaries.push_back(gradients->lower_boundaries[i].real_value - 0.5 * widths[i]);
00418 }
00419 
00420 setup(nx);
00421 
00422 if (nd > 1) {
00423 divergence.resize(nt);
00424 }
00425 }
00426 
00427 
00428 int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::real & err)
00429 {
00430 int iter = 0;
00431 
00432 if (nd == 1) {
00433 
00434 cvm::real sum = 0.0;
00435 cvm::real corr;
00436 if ( periodic[0] ) {
00437 corr = gradients->average(); // Enforce PBC by subtracting average gradient
00438 } else {
00439 corr = 0.0;
00440 }
00441 
00442 std::vector<int> ix;
00443 // Iterate over valid indices in gradient grid
00444 for (ix = new_index(); gradients->index_ok(ix); incr(ix)) {
00445 set_value(ix, sum);
00446 sum += (gradients->value_output(ix) - corr) * widths[0];
00447 }
00448 if (index_ok(ix)) {
00449 // This will happen if non-periodic: then PMF grid has one extra bin wrt gradient grid
00450 set_value(ix, sum);
00451 }
00452 
00453 } else if (nd <= 3) {
00454 
00455 nr_linbcg_sym(divergence, data, tol, itmax, iter, err);
00456 cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err) + "\n");
00457 
00458 } else {
00459 cvm::error("Cannot integrate PMF in dimension > 3\n");
00460 }
00461 
00462 return iter;
00463 }
00464 
00465 
00466 void integrate_potential::set_div()
00467 {
00468 if (nd == 1) return;
00469 for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
00470 update_div_local(ix);
00471 }
00472 }
00473 
00474 
00475 void integrate_potential::update_div_neighbors(const std::vector<int> &ix0)
00476 {
00477 std::vector<int> ix(ix0);
00478 int i, j, k;
00479 
00480 // If not periodic, expanded grid ensures that neighbors of ix0 are valid grid points
00481 if (nd == 1) {
00482 return;
00483 
00484 } else if (nd == 2) {
00485 
00486 update_div_local(ix);
00487 ix[0]++; wrap(ix);
00488 update_div_local(ix);
00489 ix[1]++; wrap(ix);
00490 update_div_local(ix);
00491 ix[0]--; wrap(ix);
00492 update_div_local(ix);
00493 
00494 } else if (nd == 3) {
00495 
00496 for (i = 0; i<2; i++) {
00497 ix[1] = ix0[1];
00498 for (j = 0; j<2; j++) {
00499 ix[2] = ix0[2];
00500 for (k = 0; k<2; k++) {
00501 wrap(ix);
00502 update_div_local(ix);
00503 ix[2]++;
00504 }
00505 ix[1]++;
00506 }
00507 ix[0]++;
00508 }
00509 }
00510 }
00511 
00512 void integrate_potential::get_grad(cvm::real * g, std::vector<int> &ix)
00513 {
00514 size_t count, i;
00515 bool edge = gradients->wrap_edge(ix); // Detect edge if non-PBC
00516 
00517 if (gradients->samples) {
00518 count = gradients->samples->value(ix);
00519 } else {
00520 count = 1;
00521 }
00522 
00523 if (!edge && count) {
00524 cvm::real const *grad = &(gradients->value(ix));
00525 cvm::real const fact = 1.0 / count;
00526 for ( i = 0; i<nd; i++ ) {
00527 g[i] = fact * grad[i];
00528 }
00529 } else {
00530 for ( i = 0; i<nd; i++ ) {
00531 g[i] = 0.0;
00532 }
00533 }
00534 }
00535 
00536 void integrate_potential::update_div_local(const std::vector<int> &ix0)
00537 {
00538 const size_t linear_index = address(ix0);
00539 int i, j, k;
00540 std::vector<int> ix = ix0;
00541 
00542 if (nd == 2) {
00543 // gradients at grid points surrounding the current scalar grid point
00544 cvm::real g00[2], g01[2], g10[2], g11[2];
00545 
00546 get_grad(g11, ix);
00547 ix[0] = ix0[0] - 1;
00548 get_grad(g01, ix);
00549 ix[1] = ix0[1] - 1;
00550 get_grad(g00, ix);
00551 ix[0] = ix0[0];
00552 get_grad(g10, ix);
00553 
00554 divergence[linear_index] = ((g10[0]-g00[0] + g11[0]-g01[0]) / widths[0]
00555 + (g01[1]-g00[1] + g11[1]-g10[1]) / widths[1]) * 0.5;
00556 } else if (nd == 3) {
00557 cvm::real gc[24]; // stores 3d gradients in 8 contiguous bins
00558 int index = 0;
00559 
00560 ix[0] = ix0[0] - 1;
00561 for (i = 0; i<2; i++) {
00562 ix[1] = ix0[1] - 1;
00563 for (j = 0; j<2; j++) {
00564 ix[2] = ix0[2] - 1;
00565 for (k = 0; k<2; k++) {
00566 get_grad(gc + index, ix);
00567 index += 3;
00568 ix[2]++;
00569 }
00570 ix[1]++;
00571 }
00572 ix[0]++;
00573 }
00574 
00575 divergence[linear_index] =
00576 ((gc[3*4]-gc[0] + gc[3*5]-gc[3*1] + gc[3*6]-gc[3*2] + gc[3*7]-gc[3*3])
00577 / widths[0]
00578 + (gc[3*2+1]-gc[0+1] + gc[3*3+1]-gc[3*1+1] + gc[3*6+1]-gc[3*4+1] + gc[3*7+1]-gc[3*5+1])
00579 / widths[1]
00580 + (gc[3*1+2]-gc[0+2] + gc[3*3+2]-gc[3*2+2] + gc[3*5+2]-gc[3*4+2] + gc[3*7+2]-gc[3*6+2])
00581 / widths[2]) * 0.25;
00582 }
00583 }
00584 
00585 
00588 void integrate_potential::atimes(const std::vector<cvm::real> &A, std::vector<cvm::real> &LA)
00589 {
00590 if (nd == 2) {
00591 // DIMENSION 2
00592 
00593 size_t index, index2;
00594 int i, j;
00595 cvm::real fact;
00596 const cvm::real ffx = 1.0 / (widths[0] * widths[0]);
00597 const cvm::real ffy = 1.0 / (widths[1] * widths[1]);
00598 const int h = nx[1];
00599 const int w = nx[0];
00600 // offsets for 4 reference points of the Laplacian stencil
00601 int xm = -h;
00602 int xp = h;
00603 int ym = -1;
00604 int yp = 1;
00605 
00606 // NOTE on performance: this version is slightly sub-optimal because
00607 // it contains two double loops on the core of the array (for x and y terms)
00608 // The slightly faster version is in commit 0254cb5a2958cb2e135f268371c4b45fad34866b
00609 // yet it is much uglier, and probably horrible to extend to dimension 3
00610 // All terms in the matrix are assigned (=) during the x loops, then updated (+=)
00611 // with the y (and z) contributions
00612 
00613 
00614 // All x components except on x edges
00615 index = h; // Skip first column
00616 
00617 // Halve the term on y edges (if any) to preserve symmetry of the Laplacian matrix
00618 // (Long Chen, Finite Difference Methods, UCI, 2017)
00619 fact = periodic[1] ? 1.0 : 0.5;
00620 
00621 for (i=1; i<w-1; i++) {
00622 // Full range of j, but factor may change on y edges (j == 0 and j == h-1)
00623 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00624 index++;
00625 for (j=1; j<h-1; j++) {
00626 LA[index] = ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00627 index++;
00628 }
00629 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00630 index++;
00631 }
00632 // Edges along x (x components only)
00633 index = 0L; // Follows left edge
00634 index2 = h * static_cast<size_t>(w - 1); // Follows right edge
00635 if (periodic[0]) {
00636 xm = h * (w - 1);
00637 xp = h;
00638 fact = periodic[1] ? 1.0 : 0.5;
00639 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00640 LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
00641 index++;
00642 index2++;
00643 for (j=1; j<h-1; j++) {
00644 LA[index] = ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00645 LA[index2] = ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
00646 index++;
00647 index2++;
00648 }
00649 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00650 LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
00651 } else {
00652 xm = -h;
00653 xp = h;
00654 fact = periodic[1] ? 1.0 : 0.5; // Halve in corners in full PBC only
00655 // lower corner, "j == 0"
00656 LA[index] = fact * ffx * (A[index + xp] - A[index]);
00657 LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
00658 index++;
00659 index2++;
00660 for (j=1; j<h-1; j++) {
00661 // x gradient (+ y term of laplacian, calculated below)
00662 LA[index] = ffx * (A[index + xp] - A[index]);
00663 LA[index2] = ffx * (A[index2 + xm] - A[index2]);
00664 index++;
00665 index2++;
00666 }
00667 // upper corner, j == h-1
00668 LA[index] = fact * ffx * (A[index + xp] - A[index]);
00669 LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
00670 }
00671 
00672 // Now adding all y components
00673 // All y components except on y edges
00674 index = 1; // Skip first element (in first row)
00675 
00676 fact = periodic[0] ? 1.0 : 0.5; // for i == 0
00677 for (i=0; i<w; i++) {
00678 // Factor of 1/2 on x edges if non-periodic
00679 if (i == 1) fact = 1.0;
00680 if (i == w - 1) fact = periodic[0] ? 1.0 : 0.5;
00681 for (j=1; j<h-1; j++) {
00682 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00683 index++;
00684 }
00685 index += 2; // skip the edges and move to next column
00686 }
00687 // Edges along y (y components only)
00688 index = 0L; // Follows bottom edge
00689 index2 = h - 1; // Follows top edge
00690 if (periodic[1]) {
00691 fact = periodic[0] ? 1.0 : 0.5;
00692 ym = h - 1;
00693 yp = 1;
00694 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00695 LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
00696 index += h;
00697 index2 += h;
00698 for (i=1; i<w-1; i++) {
00699 LA[index] += ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00700 LA[index2] += ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
00701 index += h;
00702 index2 += h;
00703 }
00704 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00705 LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
00706 } else {
00707 ym = -1;
00708 yp = 1;
00709 fact = periodic[0] ? 1.0 : 0.5; // Halve in corners in full PBC only
00710 // Left corner
00711 LA[index] += fact * ffy * (A[index + yp] - A[index]);
00712 LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
00713 index += h;
00714 index2 += h;
00715 for (i=1; i<w-1; i++) {
00716 // y gradient (+ x term of laplacian, calculated above)
00717 LA[index] += ffy * (A[index + yp] - A[index]);
00718 LA[index2] += ffy * (A[index2 + ym] - A[index2]);
00719 index += h;
00720 index2 += h;
00721 }
00722 // Right corner
00723 LA[index] += fact * ffy * (A[index + yp] - A[index]);
00724 LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
00725 }
00726 
00727 } else if (nd == 3) {
00728 // DIMENSION 3
00729 
00730 int i, j, k;
00731 size_t index, index2;
00732 cvm::real fact = 1.0;
00733 const cvm::real ffx = 1.0 / (widths[0] * widths[0]);
00734 const cvm::real ffy = 1.0 / (widths[1] * widths[1]);
00735 const cvm::real ffz = 1.0 / (widths[2] * widths[2]);
00736 const int h = nx[2]; // height
00737 const int d = nx[1]; // depth
00738 const int w = nx[0]; // width
00739 // offsets for 6 reference points of the Laplacian stencil
00740 int xm = -d * h;
00741 int xp = d * h;
00742 int ym = -h;
00743 int yp = h;
00744 int zm = -1;
00745 int zp = 1;
00746 
00747 cvm::real factx = periodic[0] ? 1 : 0.5; // factor to be applied on x edges
00748 cvm::real facty = periodic[1] ? 1 : 0.5; // same for y
00749 cvm::real factz = periodic[2] ? 1 : 0.5; // same for z
00750 cvm::real ifactx = 1 / factx;
00751 cvm::real ifacty = 1 / facty;
00752 cvm::real ifactz = 1 / factz;
00753 
00754 // All x components except on x edges
00755 index = d * static_cast<size_t>(h); // Skip left slab
00756 fact = facty * factz;
00757 for (i=1; i<w-1; i++) {
00758 for (j=0; j<d; j++) { // full range of y
00759 if (j == 1) fact *= ifacty;
00760 if (j == d-1) fact *= facty;
00761 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00762 index++;
00763 fact *= ifactz;
00764 for (k=1; k<h-1; k++) { // full range of z
00765 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00766 index++;
00767 }
00768 fact *= factz;
00769 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00770 index++;
00771 }
00772 }
00773 // Edges along x (x components only)
00774 index = 0L; // Follows left slab
00775 index2 = static_cast<size_t>(d) * h * (w - 1); // Follows right slab
00776 if (periodic[0]) {
00777 xm = d * h * (w - 1);
00778 xp = d * h;
00779 fact = facty * factz;
00780 for (j=0; j<d; j++) {
00781 if (j == 1) fact *= ifacty;
00782 if (j == d-1) fact *= facty;
00783 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00784 LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
00785 index++;
00786 index2++;
00787 fact *= ifactz;
00788 for (k=1; k<h-1; k++) {
00789 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00790 LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
00791 index++;
00792 index2++;
00793 }
00794 fact *= factz;
00795 LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
00796 LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
00797 index++;
00798 index2++;
00799 }
00800 } else {
00801 xm = -d * h;
00802 xp = d * h;
00803 fact = facty * factz;
00804 for (j=0; j<d; j++) {
00805 if (j == 1) fact *= ifacty;
00806 if (j == d-1) fact *= facty;
00807 LA[index] = fact * ffx * (A[index + xp] - A[index]);
00808 LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
00809 index++;
00810 index2++;
00811 fact *= ifactz;
00812 for (k=1; k<h-1; k++) {
00813 // x gradient (+ y, z terms of laplacian, calculated below)
00814 LA[index] = fact * ffx * (A[index + xp] - A[index]);
00815 LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
00816 index++;
00817 index2++;
00818 }
00819 fact *= factz;
00820 LA[index] = fact * ffx * (A[index + xp] - A[index]);
00821 LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
00822 index++;
00823 index2++;
00824 }
00825 }
00826 
00827 // Now adding all y components
00828 // All y components except on y edges
00829 index = h; // Skip first column (in front slab)
00830 fact = factx * factz;
00831 for (i=0; i<w; i++) { // full range of x
00832 if (i == 1) fact *= ifactx;
00833 if (i == w-1) fact *= factx;
00834 for (j=1; j<d-1; j++) {
00835 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00836 index++;
00837 fact *= ifactz;
00838 for (k=1; k<h-1; k++) {
00839 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00840 index++;
00841 }
00842 fact *= factz;
00843 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00844 index++;
00845 }
00846 index += 2 * h; // skip columns in front and back slabs
00847 }
00848 // Edges along y (y components only)
00849 index = 0L; // Follows front slab
00850 index2 = h * static_cast<size_t>(d - 1); // Follows back slab
00851 if (periodic[1]) {
00852 ym = h * (d - 1);
00853 yp = h;
00854 fact = factx * factz;
00855 for (i=0; i<w; i++) {
00856 if (i == 1) fact *= ifactx;
00857 if (i == w-1) fact *= factx;
00858 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00859 LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
00860 index++;
00861 index2++;
00862 fact *= ifactz;
00863 for (k=1; k<h-1; k++) {
00864 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00865 LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
00866 index++;
00867 index2++;
00868 }
00869 fact *= factz;
00870 LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
00871 LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
00872 index++;
00873 index2++;
00874 index += h * static_cast<size_t>(d - 1);
00875 index2 += h * static_cast<size_t>(d - 1);
00876 }
00877 } else {
00878 ym = -h;
00879 yp = h;
00880 fact = factx * factz;
00881 for (i=0; i<w; i++) {
00882 if (i == 1) fact *= ifactx;
00883 if (i == w-1) fact *= factx;
00884 LA[index] += fact * ffy * (A[index + yp] - A[index]);
00885 LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
00886 index++;
00887 index2++;
00888 fact *= ifactz;
00889 for (k=1; k<h-1; k++) {
00890 // y gradient (+ x, z terms of laplacian, calculated above and below)
00891 LA[index] += fact * ffy * (A[index + yp] - A[index]);
00892 LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
00893 index++;
00894 index2++;
00895 }
00896 fact *= factz;
00897 LA[index] += fact * ffy * (A[index + yp] - A[index]);
00898 LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
00899 index++;
00900 index2++;
00901 index += h * static_cast<size_t>(d - 1);
00902 index2 += h * static_cast<size_t>(d - 1);
00903 }
00904 }
00905 
00906 // Now adding all z components
00907 // All z components except on z edges
00908 index = 1; // Skip first element (in bottom slab)
00909 fact = factx * facty;
00910 for (i=0; i<w; i++) { // full range of x
00911 if (i == 1) fact *= ifactx;
00912 if (i == w-1) fact *= factx;
00913 for (k=1; k<h-1; k++) {
00914 LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
00915 index++;
00916 }
00917 fact *= ifacty;
00918 index += 2; // skip edge slabs
00919 for (j=1; j<d-1; j++) { // full range of y
00920 for (k=1; k<h-1; k++) {
00921 LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
00922 index++;
00923 }
00924 index += 2; // skip edge slabs
00925 }
00926 fact *= facty;
00927 for (k=1; k<h-1; k++) {
00928 LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
00929 index++;
00930 }
00931 index += 2; // skip edge slabs
00932 }
00933 // Edges along z (z components onlz)
00934 index = 0; // Follows bottom slab
00935 index2 = h - 1; // Follows top slab
00936 if (periodic[2]) {
00937 zm = h - 1;
00938 zp = 1;
00939 fact = factx * facty;
00940 for (i=0; i<w; i++) {
00941 if (i == 1) fact *= ifactx;
00942 if (i == w-1) fact *= factx;
00943 LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
00944 LA[index2] += fact * ffz * (A[index2 - zp] + A[index2 - zm] - 2.0 * A[index2]);
00945 index += h;
00946 index2 += h;
00947 fact *= ifacty;
00948 for (j=1; j<d-1; j++) {
00949 LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
00950 LA[index2] += fact * ffz * (A[index2 - zp] + A[index2 - zm] - 2.0 * A[index2]);
00951 index += h;
00952 index2 += h;
00953 }
00954 fact *= facty;
00955 LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
00956 LA[index2] += fact * ffz * (A[index2 - zp] + A[index2 - zm] - 2.0 * A[index2]);
00957 index += h;
00958 index2 += h;
00959 }
00960 } else {
00961 zm = -1;
00962 zp = 1;
00963 fact = factx * facty;
00964 for (i=0; i<w; i++) {
00965 if (i == 1) fact *= ifactx;
00966 if (i == w-1) fact *= factx;
00967 LA[index] += fact * ffz * (A[index + zp] - A[index]);
00968 LA[index2] += fact * ffz * (A[index2 + zm] - A[index2]);
00969 index += h;
00970 index2 += h;
00971 fact *= ifacty;
00972 for (j=1; j<d-1; j++) {
00973 // z gradient (+ x, y terms of laplacian, calculated above)
00974 LA[index] += fact * ffz * (A[index + zp] - A[index]);
00975 LA[index2] += fact * ffz * (A[index2 + zm] - A[index2]);
00976 index += h;
00977 index2 += h;
00978 }
00979 fact *= facty;
00980 LA[index] += fact * ffz * (A[index + zp] - A[index]);
00981 LA[index2] += fact * ffz * (A[index2 + zm] - A[index2]);
00982 index += h;
00983 index2 += h;
00984 }
00985 }
00986 }
00987 }
00988 
00989 
00990 /*
00992 void integrate_potential::asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x)
00993 {
00994 for (size_t i=0; i<int(nt); i++) {
00995 x[i] = b[i] * inv_lap_diag[i]; // Jacobi preconditioner - little benefit in tests so far
00996 }
00997 return;
00998 }*/
00999 
01000 
01001 // b : RHS of equation
01002 // x : initial guess for the solution; output is solution
01003 // itol : convergence criterion
01004 void integrate_potential::nr_linbcg_sym(const std::vector<cvm::real> &b, std::vector<cvm::real> &x, const cvm::real &tol,
01005 const int itmax, int &iter, cvm::real &err)
01006 {
01007 cvm::real ak,akden,bk,bkden,bknum,bnrm;
01008 const cvm::real EPS=1.0e-14;
01009 int j;
01010 std::vector<cvm::real> p(nt), r(nt), z(nt);
01011 
01012 iter=0;
01013 atimes(x,r);
01014 for (j=0;j<int(nt);j++) {
01015 r[j]=b[j]-r[j];
01016 }
01017 bnrm=l2norm(b);
01018 if (bnrm < EPS) {
01019 return; // Target is zero, will break relative error calc
01020 }
01021 // asolve(r,z); // precon
01022 bkden = 1.0;
01023 while (iter < itmax) {
01024 ++iter;
01025 for (bknum=0.0,j=0;j<int(nt);j++) {
01026 bknum += r[j]*r[j]; // precon: z[j]*r[j]
01027 }
01028 if (iter == 1) {
01029 for (j=0;j<int(nt);j++) {
01030 p[j] = r[j]; // precon: p[j] = z[j]
01031 }
01032 } else {
01033 bk=bknum/bkden;
01034 for (j=0;j<int(nt);j++) {
01035 p[j] = bk*p[j] + r[j]; // precon: bk*p[j] + z[j]
01036 }
01037 }
01038 bkden = bknum;
01039 atimes(p,z);
01040 for (akden=0.0,j=0;j<int(nt);j++) {
01041 akden += z[j]*p[j];
01042 }
01043 ak = bknum/akden;
01044 for (j=0;j<int(nt);j++) {
01045 x[j] += ak*p[j];
01046 r[j] -= ak*z[j];
01047 }
01048 // asolve(r,z); // precon
01049 err = l2norm(r)/bnrm;
01050 if (cvm::debug())
01051 std::cout << "iter=" << std::setw(4) << iter+1 << std::setw(12) << err << std::endl;
01052 if (err <= tol)
01053 break;
01054 }
01055 }
01056 
01057 cvm::real integrate_potential::l2norm(const std::vector<cvm::real> &x)
01058 {
01059 size_t i;
01060 cvm::real sum = 0.0;
01061 for (i=0;i<x.size();i++)
01062 sum += x[i]*x[i];
01063 return sqrt(sum);
01064 }