| version 1.16 | version 1.17 |
|---|
| |
| // -*- c++ -*- | // -*- c++ -*- |
| | |
| | // This file is part of the Collective Variables module (Colvars). |
| | // The original version of Colvars and its updates are located at: |
| | // https://github.com/colvars/colvars |
| | // Please update all Colvars source files before making any changes. |
| | // If you wish to distribute your changes, please submit them to the |
| | // Colvars repository at GitHub. |
| | |
| #ifndef COLVARGRID_H | #ifndef COLVARGRID_H |
| #define COLVARGRID_H | #define COLVARGRID_H |
| | |
| |
| inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0, | inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0, |
| int n = 0) | int n = 0) |
| { | { |
| cvm::real A0, A1; | int A0, A1, A2; |
| std::vector<int> ix; | std::vector<int> ix = ix0; |
| | |
| // factor for mesh width, 2.0 for central finite difference | |
| // but only 1.0 on edges for non-PBC coordinates | |
| cvm::real factor; | |
| | |
| if (periodic[n]) { | if (periodic[n]) { |
| factor = 2.; | |
| ix = ix0; | |
| ix[n]--; wrap(ix); | ix[n]--; wrap(ix); |
| A0 = data[address(ix)]; | A0 = data[address(ix)]; |
| ix = ix0; | ix = ix0; |
| ix[n]++; wrap(ix); | ix[n]++; wrap(ix); |
| A1 = data[address(ix)]; | A1 = data[address(ix)]; |
| | if (A0 * A1 == 0) { |
| | return 0.; // can't handle empty bins |
| } else { | } else { |
| factor = 0.; | return (std::log((cvm::real)A1) - std::log((cvm::real)A0)) |
| ix = ix0; | / (widths[n] * 2.); |
| if (ix[n] > 0) { // not left edge | |
| ix[n]--; | |
| factor += 1.; | |
| } | } |
| | } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge |
| | ix[n]--; |
| A0 = data[address(ix)]; | A0 = data[address(ix)]; |
| ix = ix0; | ix = ix0; |
| if (ix[n]+1 < nx[n]) { // not right edge | |
| ix[n]++; | ix[n]++; |
| factor += 1.; | |
| } | |
| A1 = data[address(ix)]; | A1 = data[address(ix)]; |
| } | if (A0 * A1 == 0) { |
| if (A0 == 0 || A1 == 0) { | return 0.; // can't handle empty bins |
| // can't handle empty bins | |
| return 0.; | |
| } else { | } else { |
| return (std::log((cvm::real)A1) - std::log((cvm::real)A0)) | return (std::log((cvm::real)A1) - std::log((cvm::real)A0)) |
| / (widths[n] * factor); | / (widths[n] * 2.); |
| | } |
| | } else { |
| | // edge: use 2nd order derivative |
| | int increment = (ix[n] == 0 ? 1 : -1); |
| | // move right from left edge, or the other way around |
| | A0 = data[address(ix)]; |
| | ix[n] += increment; A1 = data[address(ix)]; |
| | ix[n] += increment; A2 = data[address(ix)]; |
| | if (A0 * A1 * A2 == 0) { |
| | return 0.; // can't handle empty bins |
| | } else { |
| | return (-1.5 * std::log((cvm::real)A0) + 2. * std::log((cvm::real)A1) |
| | - 0.5 * std::log((cvm::real)A2)) * increment / widths[n]; |
| | } |
| } | } |
| } | } |
| }; | }; |