Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "HistogramBead.h"
23 : #include <vector>
24 : #include "Tools.h"
25 : #include "Keywords.h"
26 :
27 : #pragma GCC diagnostic error "-Wswitch"
28 :
29 : /*
30 : IMPORTANT NOTE FOR DEVELOPERS:
31 :
32 : If you add a new type of function in this file please add documentation for your new switching function type in function/Between.cpp
33 :
34 : */
35 : namespace PLMD {
36 :
37 : constexpr double DP2CUTOFF=6.25;
38 : //Non periodic constructor
39 719120 : HistogramBead::HistogramBead(KernelType kt,const double l,
40 : const double h,
41 719120 : const double w):
42 719120 : type(kt),
43 719120 : periodicity(Periodicity::notperiodic) {
44 719120 : set(l,h,w);
45 719120 : }
46 : //periodic constructor
47 0 : HistogramBead::HistogramBead(KernelType kt,
48 : double mlow,
49 : double mhigh,const double l,
50 : const double h,
51 0 : const double w):
52 0 : type(kt) {
53 : isPeriodic(mlow,mhigh);
54 0 : set(l,h,w);
55 0 : }
56 :
57 0 : HistogramBead::HistogramBead(const HistogramBead&)=default;
58 0 : HistogramBead::HistogramBead(HistogramBead&&)=default;
59 0 : HistogramBead& HistogramBead::operator=(const HistogramBead&)=default;
60 0 : HistogramBead& HistogramBead::operator=(HistogramBead&&)=default;
61 :
62 0 : void HistogramBead::registerKeywords( Keywords& keys ) {
63 0 : keys.add("compulsory","LOWER","the lower boundary for this particular bin");
64 0 : keys.add("compulsory","UPPER","the upper boundary for this particular bin");
65 0 : keys.add("compulsory","SMEAR","0.5","the amount to smear the Gaussian for each value in the distribution");
66 0 : }
67 :
68 58 : std::string HistogramBead::description() const {
69 58 : std::ostringstream ostr;
70 58 : switch (type) {
71 58 : case KernelType::gaussian:
72 58 : ostr<<"between "<<lowb<<" and "
73 58 : <<highb<<" width of gaussian window equals "<<width;
74 : break;
75 0 : case KernelType::triangular:
76 0 : ostr<<"between "<<lowb<<" and "
77 0 : <<highb<<" width of triangular window equals "<<width;
78 : break;
79 : }
80 58 : return ostr.str();
81 58 : }
82 :
83 0 : void HistogramBead::generateBins( const std::string& params, std::vector<std::string>& bins ) {
84 0 : std::vector<std::string> data=Tools::getWords(params);
85 0 : plumed_massert(data.size()>=1,"There is no input for this keyword");
86 :
87 0 : std::string name=data[0];
88 :
89 : unsigned nbins;
90 0 : std::vector<double> range(2);
91 : std::string smear;
92 0 : bool found_nb=Tools::parse(data,"NBINS",nbins);
93 0 : plumed_massert(found_nb,"Number of bins in histogram not found");
94 0 : bool found_r=Tools::parse(data,"LOWER",range[0]);
95 0 : plumed_massert(found_r,"Lower bound for histogram not specified");
96 0 : found_r=Tools::parse(data,"UPPER",range[1]);
97 0 : plumed_massert(found_r,"Upper bound for histogram not specified");
98 0 : plumed_massert(range[0]<range[1],"Range specification is dubious");
99 0 : bool found_b=Tools::parse(data,"SMEAR",smear);
100 0 : if(!found_b) {
101 0 : Tools::convert(0.5,smear);
102 : }
103 :
104 : std::string lb,ub;
105 0 : double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
106 0 : for(unsigned i=0; i<nbins; ++i) {
107 0 : Tools::convert( range[0]+i*delr, lb );
108 0 : Tools::convert( range[0]+(i+1)*delr, ub );
109 0 : bins.push_back( name + " " + "LOWER=" + lb + " " + "UPPER=" + ub + " " + "SMEAR=" + smear );
110 : }
111 0 : plumed_assert(bins.size()==nbins);
112 0 : }
113 :
114 58 : void HistogramBead::set( const std::string& params, std::string& errormsg ) {
115 58 : std::vector<std::string> data=Tools::getWords(params);
116 58 : if(data.size()<1) {
117 : errormsg="No input has been specified";
118 : return;
119 : }
120 :
121 58 : std::string name=data[0];
122 :
123 58 : if(name=="GAUSSIAN") {
124 58 : type=KernelType::gaussian;
125 58 : cutoff=std::sqrt(2.0*DP2CUTOFF);
126 0 : } else if(name=="TRIANGULAR") {
127 0 : type=KernelType::triangular;
128 0 : cutoff=1.0;
129 : } else {
130 0 : plumed_merror("cannot understand kernel type " + name );
131 : }
132 :
133 : double smear;
134 58 : bool found_r=Tools::parse(data,"LOWER",lowb);
135 58 : if( !found_r ) {
136 : errormsg="Lower bound has not been specified use LOWER";
137 : }
138 58 : found_r=Tools::parse(data,"UPPER",highb);
139 58 : if( !found_r ) {
140 : errormsg="Upper bound has not been specified use UPPER";
141 : }
142 58 : if( lowb>=highb ) {
143 : errormsg="Lower bound is higher than upper bound";
144 : }
145 :
146 58 : smear=0.5;
147 58 : Tools::parse(data,"SMEAR",smear);
148 58 : width=smear*(highb-lowb);
149 58 : }
150 :
151 3744806 : void HistogramBead::set(const double l, const double h, const double w) {
152 3744806 : lowb=l;
153 3744806 : highb=h;
154 3744806 : width=w;
155 3744806 : switch (type) {
156 880958 : case KernelType::gaussian : {
157 880958 : cutoff=std::sqrt(2.0*DP2CUTOFF);
158 : }
159 880958 : break;
160 2863848 : case KernelType::triangular : {
161 2863848 : cutoff=1.0;
162 : }
163 : }
164 3744806 : }
165 :
166 78 : HistogramBead::KernelType HistogramBead::getKernelType( const std::string& ktype ) {
167 78 : if(ktype=="gaussian") {
168 : return KernelType::gaussian;
169 5 : } else if(ktype=="triangular") {
170 : return KernelType::triangular;
171 : } else {
172 0 : plumed_merror("cannot understand kernel type " + ktype );
173 : }
174 : }
175 :
176 0 : void HistogramBead::setKernelType( const std::string& ktype ) {
177 0 : if(ktype=="gaussian") {
178 0 : type=KernelType::gaussian;
179 0 : } else if(ktype=="triangular") {
180 0 : type=KernelType::triangular;
181 : } else {
182 0 : plumed_merror("cannot understand kernel type " + ktype );
183 : }
184 0 : }
185 :
186 0 : void HistogramBead::setKernelType( KernelType ktype ) {
187 0 : type=ktype;
188 0 : }
189 :
190 397966 : double HistogramBead::calculate(const double x, double& df ) const {
191 : double res=0.0;
192 397966 : switch (type) {
193 397966 : case KernelType::gaussian : {
194 397966 : const double lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
195 397966 : const double upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
196 397966 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(PLMD::twopi)*width );
197 397966 : res = 0.5*( erf( upperB ) - erf( lowB ) );
198 : }
199 397966 : break;
200 0 : case KernelType::triangular : {
201 0 : const double lowB = ( difference( x, lowb ) / width );
202 0 : const double upperB = ( difference( x, highb ) / width );
203 0 : df=0.0;
204 0 : if( std::fabs(lowB)<1. ) {
205 0 : df = (1 - std::fabs(lowB)) / width;
206 : }
207 0 : if( std::fabs(upperB)<1. ) {
208 0 : df -= (1 - std::fabs(upperB)) / width;
209 : }
210 0 : if (upperB<=-1. || lowB >=1.) {
211 : break;
212 : }
213 0 : const double ia = (lowB>-1.0) ? lowB : -1.0;
214 0 : const double ib = (upperB<1.0) ? upperB : 1.0;
215 0 : res = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
216 : }
217 : }
218 397966 : return res;
219 : }
220 :
221 425128 : double HistogramBead::calculateWithCutoff( double x, double& df ) const {
222 425128 : return calculateWithCutoff( x, lowb, highb, width, df );
223 : }
224 :
225 958728 : double HistogramBead::calculateWithCutoff(double x,
226 : double l,
227 : double u,
228 : double s,
229 : double& df ) const {
230 958728 : double lowB = difference( x, l ) / s;
231 958728 : double upperB = difference( x, u ) / s;
232 958728 : if( upperB<=-cutoff || lowB>=cutoff ) {
233 87679 : df=0.0;
234 87679 : return 0.0;
235 : }
236 : double f=0.0;
237 871049 : switch (type) {
238 446863 : case KernelType::gaussian : {
239 446863 : lowB /= std::sqrt(2.0);
240 446863 : upperB /= std::sqrt(2.0);
241 446863 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(2*pi)*s );
242 446863 : f = 0.5*( erf( upperB ) - erf( lowB ) );
243 : }
244 446863 : break;
245 424186 : case KernelType::triangular : {
246 424186 : df=0;
247 424186 : if( std::fabs(lowB)<1. ) {
248 108180 : df = (1.0 - std::fabs(lowB)) / s;
249 : }
250 424186 : if( std::fabs(upperB)<1. ) {
251 108180 : df -= (1.0 - std::fabs(upperB)) / s;
252 : }
253 424186 : if (upperB<=-1. || lowB >=1.) {
254 : f=0.0;
255 : } else {
256 : double ia, ib;
257 424186 : if( lowB>-1.0 ) {
258 : ia=lowB;
259 : } else {
260 : ia=-1.0;
261 : }
262 424186 : if( upperB<1.0 ) {
263 : ib=upperB;
264 : } else {
265 : ib=1.0;
266 : }
267 424186 : f = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
268 : }
269 : }
270 : }
271 : return f;
272 : }
273 :
274 9360 : double HistogramBead::lboundDerivative( const double x ) const {
275 9360 : switch (type) {
276 9360 : case KernelType::gaussian : {
277 9360 : double lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
278 9360 : return exp( -lowB*lowB ) / ( std::sqrt(2*pi)*width );
279 : }
280 : break;
281 : case KernelType::triangular : {
282 : assert(false);
283 : //TODO: something about this:
284 : // double lowB = fabs( difference( x, lowb ) / width );
285 : // if( lowB<1 ) {
286 : // return ( 1 - (lowB) ) / 2*width;
287 : // } else {
288 : // return 0;
289 : // }
290 : }
291 : }
292 : //warning::unreachable
293 : return 0.0;
294 : }
295 :
296 9360 : double HistogramBead::uboundDerivative( const double x ) const {
297 9360 : switch (type) {
298 9360 : case KernelType::gaussian : {
299 9360 : double upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
300 9360 : return exp( -upperB*upperB ) / ( std::sqrt(2*pi)*width );
301 : }
302 : break;
303 : case KernelType::triangular : {
304 : assert(false);
305 : //TODO: something about this:
306 : // double upperB = fabs( difference( x, highb ) / width );
307 : // if( upperB<1 ) {
308 : // return ( 1 - (upperB) ) / 2*width;
309 : // } else {
310 : // return 0;
311 : // }
312 : }
313 : }
314 : //warning::unreachable
315 : return 0.0;
316 : }
317 :
318 : inline
319 2732108 : double HistogramBead::difference( double d1, const double d2 ) const {
320 2732108 : switch (periodicity) {
321 : case Periodicity::notperiodic :
322 : break;
323 1067242 : case Periodicity::periodic: {
324 : // Make sure the point is in the target range
325 1067242 : d1*=inv_max_minus_min;
326 1067242 : d1=Tools::pbc(d1);
327 1067242 : d1*=max_minus_min;
328 : }
329 : }
330 2732108 : return d2-d1;
331 : }
332 :
333 : }
|
