Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2021 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 :
39 : $$
40 : p(\mathbf{s}) = \sum_{i} \, w_{i}
41 : \prod_{k}^{d}
42 : \frac{\lambda_{k,i}}{2}
43 : \,
44 : \exp\left[
45 : \frac{\lambda_{k,i}}{2}
46 : (2 \mu_{k,i} + \lambda_{k,i} \sigma_{k,i}^2 -2 s_{k})
47 : \right]
48 : \,
49 : \mathrm{erfc}\left[
50 : \frac{\mu_{k,i} + \lambda_{k,i} \sigma_{k,i}^2 - s_{k})}{\sqrt{2} \sigma_{k,i}}
51 : \right]
52 : $$
53 :
54 : where $(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})$
55 : are the centers of the Gaussian component,
56 : $(\sigma_{1,i},\sigma_{2,i},\ldots,\sigma_{d,i})$ are the
57 : standard deviations of the Gaussian component,
58 : $(\lambda_{1,i},\lambda_{2,i},\ldots,\lambda_{d,i})$ are the
59 : rate parameters of the exponential component, and
60 : $\mathrm{erfc}(x)=1-\mathrm{erf}(x)$ is the
61 : complementary error function.
62 : The weights $w_{i}$ are normalized to 1, $\sum_{i}w_{i}=1$.
63 :
64 : The centers $(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})$ are
65 : given using the numbered CENTER keywords, the standard deviations
66 : $(\sigma_{1,i},\sigma_{2,i},\ldots,\sigma_{d,i})$ using the
67 : the numbered SIGMA keywords, and the rate parameters
68 : $(\lambda_{1,i},\lambda_{2,i},\ldots,\lambda_{d,i})$ using the
69 : numbered LAMBDA keywords.
70 : The weights are given using the WEIGHTS keywords, if no weights are
71 : given are all terms weighted equally.
72 :
73 : ## Examples
74 :
75 : An exponentially modified Gaussian distribution in one-dimension
76 :
77 : ```plumed
78 : td1: TD_EXPONENTIALLY_MODIFIED_GAUSSIAN CENTER1=-10.0 SIGMA1=1.0 LAMBDA1=0.25
79 : ```
80 :
81 : A sum of two one-dimensional exponentially modified Gaussian distributions
82 :
83 : ```plumed
84 : TD_EXPONENTIALLY_MODIFIED_GAUSSIAN ...
85 : CENTER1=-10.0 SIGMA1=1.0 LAMBDA1=0.5
86 : CENTER2=+10.0 SIGMA2=1.0 LAMBDA2=1.0
87 : WEIGHTS=2.0,1.0
88 : LABEL=td1
89 : ... TD_EXPONENTIALLY_MODIFIED_GAUSSIAN
90 : ```
91 :
92 : A sum of two two-dimensional exponentially modified Gaussian distributions
93 :
94 : ```plumed
95 : TD_EXPONENTIALLY_MODIFIED_GAUSSIAN ...
96 : CENTER1=-5.0,+5.0 SIGMA1=1.0,1.0 LAMBDA1=0.5,0.5
97 : CENTER2=+5.0,+5.0 SIGMA2=1.0,1.0 LAMBDA2=1.0,1.0
98 : WEIGHTS=1.0,1.0
99 : LABEL=td1
100 : ... TD_EXPONENTIALLY_MODIFIED_GAUSSIAN
101 : ```
102 :
103 :
104 :
105 :
106 :
107 : */
108 : //+ENDPLUMEDOC
109 :
110 : class TD_ExponentiallyModifiedGaussian: public TargetDistribution {
111 : std::vector< std::vector<double> > centers_;
112 : std::vector< std::vector<double> > sigmas_;
113 : std::vector< std::vector<double> > lambdas_;
114 : std::vector<double> weights_;
115 : unsigned int ncenters_;
116 : double ExponentiallyModifiedGaussianDiagonal(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&) const;
117 : public:
118 : static void registerKeywords(Keywords&);
119 : explicit TD_ExponentiallyModifiedGaussian(const ActionOptions& ao);
120 : double getValue(const std::vector<double>&) const override;
121 : };
122 :
123 :
124 : PLUMED_REGISTER_ACTION(TD_ExponentiallyModifiedGaussian,"TD_EXPONENTIALLY_MODIFIED_GAUSSIAN")
125 :
126 :
127 8 : void TD_ExponentiallyModifiedGaussian::registerKeywords(Keywords& keys) {
128 8 : TargetDistribution::registerKeywords(keys);
129 8 : keys.add("numbered","CENTER","The center of each exponentially modified Gaussian distributions.");
130 8 : keys.add("numbered","SIGMA","The sigma parameters for each exponentially modified Gaussian distributions.");
131 8 : keys.add("numbered","LAMBDA","The lambda parameters for each exponentially modified Gaussian distributions");
132 8 : keys.add("optional","WEIGHTS","The weights of the distributions. By default all are weighted equally.");
133 8 : keys.use("WELLTEMPERED_FACTOR");
134 8 : keys.use("SHIFT_TO_ZERO");
135 8 : keys.use("NORMALIZE");
136 8 : }
137 :
138 :
139 6 : TD_ExponentiallyModifiedGaussian::TD_ExponentiallyModifiedGaussian(const ActionOptions& ao):
140 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
141 12 : centers_(0),
142 6 : sigmas_(0),
143 6 : lambdas_(0),
144 6 : weights_(0),
145 12 : ncenters_(0) {
146 9 : for(unsigned int i=1;; i++) {
147 : std::vector<double> tmp_center;
148 30 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {
149 : break;
150 : }
151 9 : centers_.push_back(tmp_center);
152 9 : }
153 9 : for(unsigned int i=1;; i++) {
154 : std::vector<double> tmp_sigma;
155 30 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {
156 : break;
157 : }
158 20 : for(unsigned int k=0; k<tmp_sigma.size(); k++) {
159 11 : if(tmp_sigma[k]<=0.0) {
160 0 : plumed_merror(getName()+": the values given in SIGMA should be positive");
161 : }
162 : }
163 9 : sigmas_.push_back(tmp_sigma);
164 9 : }
165 9 : for(unsigned int i=1;; i++) {
166 : std::vector<double> tmp_lambda;
167 30 : if(!parseNumberedVector("LAMBDA",i,tmp_lambda) ) {
168 : break;
169 : }
170 20 : for(unsigned int k=0; k<tmp_lambda.size(); k++) {
171 11 : if(tmp_lambda[k]<=0.0) {
172 0 : plumed_merror(getName()+": the values given in LAMBDA should be positive");
173 : }
174 : }
175 9 : lambdas_.push_back(tmp_lambda);
176 9 : }
177 : //
178 6 : if(centers_.size()==0) {
179 0 : plumed_merror(getName()+": CENTER keywords seem to be missing. Note that numbered keywords start at CENTER1.");
180 : }
181 : //
182 6 : if(centers_.size()!=sigmas_.size() || centers_.size()!=lambdas_.size() ) {
183 0 : plumed_merror(getName()+": there has to be an equal amount of CENTER, SIGMA, and LAMBDA keywords");
184 : }
185 : //
186 6 : setDimension(centers_[0].size());
187 6 : ncenters_ = centers_.size();
188 : //
189 : // check centers and sigmas
190 15 : for(unsigned int i=0; i<ncenters_; i++) {
191 9 : if(centers_[i].size()!=getDimension()) {
192 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
193 : }
194 9 : if(sigmas_[i].size()!=getDimension()) {
195 0 : plumed_merror(getName()+": one of the SIGMA keyword does not match the given dimension");
196 : }
197 9 : if(lambdas_[i].size()!=getDimension()) {
198 0 : plumed_merror(getName()+": one of the LAMBDA keyword does not match the given dimension");
199 : }
200 : }
201 : //
202 12 : parseVector("WEIGHTS",weights_);
203 6 : if(weights_.size()==0) {
204 4 : weights_.assign(centers_.size(),1.0);
205 : }
206 6 : if(centers_.size()!=weights_.size()) {
207 0 : plumed_merror(getName()+": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
208 : }
209 : //
210 : double sum_weights=0.0;
211 15 : for(unsigned int i=0; i<weights_.size(); i++) {
212 9 : sum_weights+=weights_[i];
213 : }
214 15 : for(unsigned int i=0; i<weights_.size(); i++) {
215 9 : weights_[i]/=sum_weights;
216 : }
217 : //
218 6 : checkRead();
219 6 : }
220 :
221 :
222 11206 : double TD_ExponentiallyModifiedGaussian::getValue(const std::vector<double>& argument) const {
223 : double value=0.0;
224 33015 : for(unsigned int i=0; i<ncenters_; i++) {
225 21809 : value+=weights_[i]*ExponentiallyModifiedGaussianDiagonal(argument,centers_[i],sigmas_[i],lambdas_[i]);
226 : }
227 11206 : return value;
228 : }
229 :
230 :
231 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 {
232 : double value = 1.0;
233 64020 : for(unsigned int k=0; k<argument.size(); k++) {
234 42211 : double arg1 = 0.5*lambda[k]*(2.0*center[k]+lambda[k]*sigma[k]*sigma[k]-2.0*argument[k]);
235 42211 : double arg2 = (center[k]+lambda[k]*sigma[k]*sigma[k]-argument[k])/(sqrt(2.0)*sigma[k]);
236 42211 : value *= 0.5*lambda[k]*exp(arg1)*erfc(arg2);
237 : }
238 21809 : return value;
239 : }
240 :
241 :
242 :
243 : }
244 : }
|
