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 <limits>
25 : #include "Tools.h"
26 : #include "Keywords.h"
27 :
28 : /*
29 : IMPORTANT NOTE FOR DEVELOPERS:
30 :
31 : If you add a new type of function in this file please add documentation for your new switching function type in function/Between.cpp
32 :
33 : */
34 :
35 : namespace PLMD {
36 :
37 0 : void HistogramBead::registerKeywords( Keywords& keys ) {
38 0 : keys.add("compulsory","LOWER","the lower boundary for this particular bin");
39 0 : keys.add("compulsory","UPPER","the upper boundary for this particular bin");
40 0 : keys.add("compulsory","SMEAR","0.5","the amount to smear the Gaussian for each value in the distribution");
41 0 : }
42 :
43 260469 : HistogramBead::HistogramBead():
44 260469 : init(false),
45 260469 : lowb(0.0),
46 260469 : highb(0.0),
47 260469 : width(0.0),
48 260469 : cutoff(std::numeric_limits<double>::max()),
49 260469 : type(gaussian),
50 260469 : periodicity(unset),
51 260469 : min(0.0),
52 260469 : max(0.0),
53 260469 : max_minus_min(0.0),
54 260469 : inv_max_minus_min(0.0) {
55 260469 : }
56 :
57 53 : std::string HistogramBead::description() const {
58 53 : std::ostringstream ostr;
59 53 : ostr<<"between "<<lowb<<" and "<<highb<<" width of gaussian window equals "<<width;
60 53 : return ostr.str();
61 53 : }
62 :
63 0 : void HistogramBead::generateBins( const std::string& params, std::vector<std::string>& bins ) {
64 0 : std::vector<std::string> data=Tools::getWords(params);
65 0 : plumed_massert(data.size()>=1,"There is no input for this keyword");
66 :
67 0 : std::string name=data[0];
68 :
69 : unsigned nbins;
70 0 : std::vector<double> range(2);
71 : std::string smear;
72 0 : bool found_nb=Tools::parse(data,"NBINS",nbins);
73 0 : plumed_massert(found_nb,"Number of bins in histogram not found");
74 0 : bool found_r=Tools::parse(data,"LOWER",range[0]);
75 0 : plumed_massert(found_r,"Lower bound for histogram not specified");
76 0 : found_r=Tools::parse(data,"UPPER",range[1]);
77 0 : plumed_massert(found_r,"Upper bound for histogram not specified");
78 0 : plumed_massert(range[0]<range[1],"Range specification is dubious");
79 0 : bool found_b=Tools::parse(data,"SMEAR",smear);
80 0 : if(!found_b) {
81 0 : Tools::convert(0.5,smear);
82 : }
83 :
84 : std::string lb,ub;
85 0 : double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
86 0 : for(unsigned i=0; i<nbins; ++i) {
87 0 : Tools::convert( range[0]+i*delr, lb );
88 0 : Tools::convert( range[0]+(i+1)*delr, ub );
89 0 : bins.push_back( name + " " + "LOWER=" + lb + " " + "UPPER=" + ub + " " + "SMEAR=" + smear );
90 : }
91 0 : plumed_assert(bins.size()==nbins);
92 0 : }
93 :
94 53 : void HistogramBead::set( const std::string& params, std::string& errormsg ) {
95 53 : std::vector<std::string> data=Tools::getWords(params);
96 53 : if(data.size()<1) {
97 : errormsg="No input has been specified";
98 : return;
99 : }
100 :
101 53 : std::string name=data[0];
102 : const double DP2CUTOFF=6.25;
103 53 : if(name=="GAUSSIAN") {
104 53 : type=gaussian;
105 53 : cutoff=std::sqrt(2.0*DP2CUTOFF);
106 0 : } else if(name=="TRIANGULAR") {
107 0 : type=triangular;
108 0 : cutoff=1.;
109 : } else {
110 0 : plumed_merror("cannot understand kernel type " + name );
111 : }
112 :
113 : double smear;
114 53 : bool found_r=Tools::parse(data,"LOWER",lowb);
115 53 : if( !found_r ) {
116 : errormsg="Lower bound has not been specified use LOWER";
117 : }
118 53 : found_r=Tools::parse(data,"UPPER",highb);
119 53 : if( !found_r ) {
120 : errormsg="Upper bound has not been specified use UPPER";
121 : }
122 53 : if( lowb>=highb ) {
123 : errormsg="Lower bound is higher than upper bound";
124 : }
125 :
126 53 : smear=0.5;
127 53 : Tools::parse(data,"SMEAR",smear);
128 53 : width=smear*(highb-lowb);
129 53 : init=true;
130 53 : }
131 :
132 2668474 : void HistogramBead::set( double l, double h, double w) {
133 2668474 : init=true;
134 2668474 : lowb=l;
135 2668474 : highb=h;
136 2668474 : width=w;
137 : const double DP2CUTOFF=6.25;
138 2668474 : if( type==gaussian ) {
139 733122 : cutoff=std::sqrt(2.0*DP2CUTOFF);
140 1935352 : } else if( type==triangular ) {
141 1935352 : cutoff=1.;
142 : } else {
143 0 : plumed_error();
144 : }
145 2668474 : }
146 :
147 260148 : void HistogramBead::setKernelType( const std::string& ktype ) {
148 260148 : if(ktype=="gaussian") {
149 233896 : type=gaussian;
150 26252 : } else if(ktype=="triangular") {
151 26252 : type=triangular;
152 : } else {
153 0 : plumed_merror("cannot understand kernel type " + ktype );
154 : }
155 260148 : }
156 :
157 251184 : double HistogramBead::calculate( double x, double& df ) const {
158 : plumed_dbg_assert(init && periodicity!=unset );
159 : double lowB, upperB, f;
160 251184 : if( type==gaussian ) {
161 251184 : lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
162 251184 : upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
163 251184 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(2*pi)*width );
164 251184 : f = 0.5*( erf( upperB ) - erf( lowB ) );
165 0 : } else if( type==triangular ) {
166 0 : lowB = ( difference( x, lowb ) / width );
167 0 : upperB = ( difference( x, highb ) / width );
168 0 : df=0;
169 0 : if( std::fabs(lowB)<1. ) {
170 0 : df = (1 - std::fabs(lowB)) / width;
171 : }
172 0 : if( std::fabs(upperB)<1. ) {
173 0 : df -= (1 - std::fabs(upperB)) / width;
174 : }
175 0 : if (upperB<=-1. || lowB >=1.) {
176 : f=0.;
177 : } else {
178 : double ia, ib;
179 0 : if( lowB>-1.0 ) {
180 : ia=lowB;
181 : } else {
182 : ia=-1.0;
183 : }
184 0 : if( upperB<1.0 ) {
185 : ib=upperB;
186 : } else {
187 : ib=1.0;
188 : }
189 0 : f = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
190 : }
191 : } else {
192 0 : plumed_merror("function type does not exist");
193 : }
194 251184 : return f;
195 : }
196 :
197 920596 : double HistogramBead::calculateWithCutoff( double x, double& df ) const {
198 : plumed_dbg_assert(init && periodicity!=unset );
199 :
200 : double lowB, upperB, f;
201 920596 : lowB = difference( x, lowb ) / width ;
202 920596 : upperB = difference( x, highb ) / width;
203 920596 : if( upperB<=-cutoff || lowB>=cutoff ) {
204 87679 : df=0;
205 87679 : return 0;
206 : }
207 :
208 832917 : if( type==gaussian ) {
209 446395 : lowB /= std::sqrt(2.0);
210 446395 : upperB /= std::sqrt(2.0);
211 446395 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(2*pi)*width );
212 446395 : f = 0.5*( erf( upperB ) - erf( lowB ) );
213 386522 : } else if( type==triangular ) {
214 386522 : df=0;
215 386522 : if( std::fabs(lowB)<1. ) {
216 98488 : df = (1 - std::fabs(lowB)) / width;
217 : }
218 386522 : if( std::fabs(upperB)<1. ) {
219 98488 : df -= (1 - std::fabs(upperB)) / width;
220 : }
221 386522 : if (upperB<=-1. || lowB >=1.) {
222 : f=0.;
223 : } else {
224 : double ia, ib;
225 386522 : if( lowB>-1.0 ) {
226 : ia=lowB;
227 : } else {
228 : ia=-1.0;
229 : }
230 386522 : if( upperB<1.0 ) {
231 : ib=upperB;
232 : } else {
233 : ib=1.0;
234 : }
235 386522 : f = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
236 : }
237 : } else {
238 0 : plumed_merror("function type does not exist");
239 : }
240 : return f;
241 : }
242 :
243 4680 : double HistogramBead::lboundDerivative( const double& x ) const {
244 4680 : if( type==gaussian ) {
245 4680 : double lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
246 4680 : return exp( -lowB*lowB ) / ( std::sqrt(2*pi)*width );
247 0 : } else if ( type==triangular ) {
248 0 : plumed_error();
249 : // lowB = fabs( difference( x, lowb ) / width );
250 : // if( lowB<1 ) return ( 1 - (lowB) ) / 2*width;
251 : // else return 0;
252 : } else {
253 0 : plumed_merror("function type does not exist");
254 : }
255 : return 0;
256 : }
257 :
258 4680 : double HistogramBead::uboundDerivative( const double& x ) const {
259 : plumed_dbg_assert(init && periodicity!=unset );
260 4680 : if( type==gaussian ) {
261 4680 : double upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
262 4680 : return exp( -upperB*upperB ) / ( std::sqrt(2*pi)*width );
263 0 : } else if ( type==triangular ) {
264 0 : plumed_error();
265 : // upperB = fabs( difference( x, highb ) / width );
266 : // if( upperB<1 ) return ( 1 - (upperB) ) / 2*width;
267 : // else return 0;
268 : } else {
269 0 : plumed_merror("function type does not exist");
270 : }
271 : return 0;
272 : }
273 :
274 : }
|
