Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module 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 : The VES code module 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 the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "TargetDistribution.h"
24 :
25 : #include "core/ActionRegister.h"
26 :
27 :
28 : namespace PLMD {
29 : namespace ves {
30 :
31 : //+PLUMEDOC VES_TARGETDIST TD_EXPONENTIALLY_MODIFIED_GAUSSIAN
32 : /*
33 : Target distribution given by a sum of exponentially modified Gaussian distributions (static).
34 :
35 : Employ a target distribution that is given by a sum where each
36 : term is a product of one-dimensional
37 : [exponentially modified Gaussian distributions](http://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution),
38 : \f[
39 : p(\mathbf{s}) = \sum_{i} \, w_{i}
40 : \prod_{k}^{d}
41 : \frac{\lambda_{k,i}}{2}
42 : \,
43 : \exp\left[
44 : \frac{\lambda_{k,i}}{2}
45 : (2 \mu_{k,i} + \lambda_{k,i} \sigma_{k,i}^2 -2 s_{k})
46 : \right]
47 : \,
48 : \mathrm{erfc}\left[
49 : \frac{\mu_{k,i} + \lambda_{k,i} \sigma_{k,i}^2 - s_{k})}{\sqrt{2} \sigma_{k,i}}
50 : \right]
51 : \f]
52 : where \f$(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})\f$
53 : are the centers of the Gaussian component,
54 : \f$(\sigma_{1,i},\sigma_{2,i},\ldots,\sigma_{d,i})\f$ are the
55 : standard deviations of the Gaussian component,
56 : \f$(\lambda_{1,i},\lambda_{2,i},\ldots,\lambda_{d,i})\f$ are the
57 : rate parameters of the exponential component, and
58 : \f$\mathrm{erfc}(x)=1-\mathrm{erf}(x)\f$ is the
59 : complementary error function.
60 : The weights \f$w_{i}\f$ are normalized to 1, \f$\sum_{i}w_{i}=1\f$.
61 :
62 : The centers \f$(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})\f$ are
63 : given using the numbered CENTER keywords, the standard deviations
64 : \f$(\sigma_{1,i},\sigma_{2,i},\ldots,\sigma_{d,i})\f$ using the
65 : the numbered SIGMA keywords, and the rate parameters
66 : \f$(\lambda_{1,i},\lambda_{2,i},\ldots,\lambda_{d,i})\f$ using the
67 : numbered LAMBDA keywords.
68 : The weights are given using the WEIGHTS keywords, if no weights are
69 : given are all terms weighted equally.
70 :
71 : \par Examples
72 :
73 : An exponentially modified Gaussian distribution in one-dimension
74 : \plumedfile
75 : td1: TD_EXPONENTIALLY_MODIFIED_GAUSSIAN CENTER1=-10.0 SIGMA1=1.0 LAMBDA1=0.25
76 : \endplumedfile
77 :
78 : A sum of two one-dimensional exponentially modified Gaussian distributions
79 : \plumedfile
80 : TD_EXPONENTIALLY_MODIFIED_GAUSSIAN ...
81 : CENTER1=-10.0 SIGMA1=1.0 LAMBDA1=0.5
82 : CENTER2=+10.0 SIGMA2=1.0 LAMBDA2=1.0
83 : WEIGHTS=2.0,1.0
84 : LABEL=td1
85 : ... TD_EXPONENTIALLY_MODIFIED_GAUSSIAN
86 : \endplumedfile
87 :
88 : A sum of two two-dimensional exponentially modified Gaussian distributions
89 : \plumedfile
90 : TD_EXPONENTIALLY_MODIFIED_GAUSSIAN ...
91 : CENTER1=-5.0,+5.0 SIGMA1=1.0,1.0 LAMBDA1=0.5,0.5
92 : CENTER2=+5.0,+5.0 SIGMA2=1.0,1.0 LAMBDA2=1.0,1.0
93 : WEIGHTS=1.0,1.0
94 : LABEL=td1
95 : ... TD_EXPONENTIALLY_MODIFIED_GAUSSIAN
96 : \endplumedfile
97 :
98 :
99 :
100 :
101 :
102 : */
103 : //+ENDPLUMEDOC
104 :
105 18 : class TD_ExponentiallyModifiedGaussian: public TargetDistribution {
106 : std::vector< std::vector<double> > centers_;
107 : std::vector< std::vector<double> > sigmas_;
108 : std::vector< std::vector<double> > lambdas_;
109 : std::vector<double> weights_;
110 : unsigned int ncenters_;
111 : double ExponentiallyModifiedGaussianDiagonal(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&) const;
112 : public:
113 : static void registerKeywords(Keywords&);
114 : explicit TD_ExponentiallyModifiedGaussian(const ActionOptions& ao);
115 : double getValue(const std::vector<double>&) const;
116 : };
117 :
118 :
119 6458 : PLUMED_REGISTER_ACTION(TD_ExponentiallyModifiedGaussian,"TD_EXPONENTIALLY_MODIFIED_GAUSSIAN")
120 :
121 :
122 7 : void TD_ExponentiallyModifiedGaussian::registerKeywords(Keywords& keys) {
123 7 : TargetDistribution::registerKeywords(keys);
124 28 : keys.add("numbered","CENTER","The center of each exponentially modified Gaussian distributions.");
125 28 : keys.add("numbered","SIGMA","The sigma parameters for each exponentially modified Gaussian distributions.");
126 28 : keys.add("numbered","LAMBDA","The lambda parameters for each exponentially modified Gaussian distributions");
127 28 : keys.add("optional","WEIGHTS","The weights of the distributions. By default all are weighted equally.");
128 14 : keys.use("WELLTEMPERED_FACTOR");
129 14 : keys.use("SHIFT_TO_ZERO");
130 14 : keys.use("NORMALIZE");
131 7 : }
132 :
133 :
134 6 : TD_ExponentiallyModifiedGaussian::TD_ExponentiallyModifiedGaussian(const ActionOptions& ao):
135 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
136 : centers_(0),
137 : sigmas_(0),
138 : lambdas_(0),
139 : weights_(0),
140 6 : ncenters_(0)
141 : {
142 9 : for(unsigned int i=1;; i++) {
143 : std::vector<double> tmp_center;
144 30 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {break;}
145 9 : centers_.push_back(tmp_center);
146 9 : }
147 9 : for(unsigned int i=1;; i++) {
148 : std::vector<double> tmp_sigma;
149 30 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {break;}
150 51 : for(unsigned int k=0; k<tmp_sigma.size(); k++) {
151 11 : if(tmp_sigma[k]<=0.0) {plumed_merror(getName()+": the values given in SIGMA should be postive");}
152 : }
153 9 : sigmas_.push_back(tmp_sigma);
154 9 : }
155 9 : for(unsigned int i=1;; i++) {
156 : std::vector<double> tmp_lambda;
157 30 : if(!parseNumberedVector("LAMBDA",i,tmp_lambda) ) {break;}
158 51 : for(unsigned int k=0; k<tmp_lambda.size(); k++) {
159 11 : if(tmp_lambda[k]<=0.0) {plumed_merror(getName()+": the values given in LAMBDA should be postive");}
160 : }
161 9 : lambdas_.push_back(tmp_lambda);
162 9 : }
163 : //
164 6 : if(centers_.size()==0) {
165 0 : plumed_merror(getName()+": CENTER keywords seem to be missing. Note that numbered keywords start at CENTER1.");
166 : }
167 : //
168 12 : if(centers_.size()!=sigmas_.size() || centers_.size()!=lambdas_.size() ) {
169 0 : plumed_merror(getName()+": there has to be an equal amount of CENTER, SIGMA, and LAMBDA keywords");
170 : }
171 : //
172 6 : setDimension(centers_[0].size());
173 6 : ncenters_ = centers_.size();
174 : //
175 : // check centers and sigmas
176 24 : for(unsigned int i=0; i<ncenters_; i++) {
177 18 : if(centers_[i].size()!=getDimension()) {
178 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
179 : }
180 9 : if(sigmas_[i].size()!=getDimension()) {
181 0 : plumed_merror(getName()+": one of the SIGMA keyword does not match the given dimension");
182 : }
183 9 : if(lambdas_[i].size()!=getDimension()) {
184 0 : plumed_merror(getName()+": one of the LAMBDA keyword does not match the given dimension");
185 : }
186 : }
187 : //
188 12 : parseVector("WEIGHTS",weights_);
189 10 : if(weights_.size()==0) {weights_.assign(centers_.size(),1.0);}
190 6 : if(centers_.size()!=weights_.size()) {
191 0 : plumed_merror(getName()+": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
192 : }
193 : //
194 : double sum_weights=0.0;
195 24 : for(unsigned int i=0; i<weights_.size(); i++) {sum_weights+=weights_[i];}
196 39 : for(unsigned int i=0; i<weights_.size(); i++) {weights_[i]/=sum_weights;}
197 : //
198 6 : checkRead();
199 6 : }
200 :
201 :
202 11206 : double TD_ExponentiallyModifiedGaussian::getValue(const std::vector<double>& argument) const {
203 : double value=0.0;
204 54824 : for(unsigned int i=0; i<ncenters_; i++) {
205 65427 : value+=weights_[i]*ExponentiallyModifiedGaussianDiagonal(argument,centers_[i],sigmas_[i],lambdas_[i]);
206 : }
207 11206 : return value;
208 : }
209 :
210 :
211 21809 : double TD_ExponentiallyModifiedGaussian::ExponentiallyModifiedGaussianDiagonal(const std::vector<double>& argument, const std::vector<double>& center, const std::vector<double>& sigma, const std::vector<double>& lambda) const {
212 : double value = 1.0;
213 170251 : for(unsigned int k=0; k<argument.size(); k++) {
214 168844 : double arg1 = 0.5*lambda[k]*(2.0*center[k]+lambda[k]*sigma[k]*sigma[k]-2.0*argument[k]);
215 42211 : double arg2 = (center[k]+lambda[k]*sigma[k]*sigma[k]-argument[k])/(sqrt(2.0)*sigma[k]);
216 42211 : value *= 0.5*lambda[k]*exp(arg1)*erfc(arg2);
217 : }
218 21809 : return value;
219 : }
220 :
221 :
222 :
223 : }
224 4839 : }
|
