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 : \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 : 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 override;
116 : };
117 :
118 :
119 : PLUMED_REGISTER_ACTION(TD_ExponentiallyModifiedGaussian,"TD_EXPONENTIALLY_MODIFIED_GAUSSIAN")
120 :
121 :
122 8 : void TD_ExponentiallyModifiedGaussian::registerKeywords(Keywords& keys) {
123 8 : TargetDistribution::registerKeywords(keys);
124 8 : keys.add("numbered","CENTER","The center of each exponentially modified Gaussian distributions.");
125 8 : keys.add("numbered","SIGMA","The sigma parameters for each exponentially modified Gaussian distributions.");
126 8 : keys.add("numbered","LAMBDA","The lambda parameters for each exponentially modified Gaussian distributions");
127 8 : keys.add("optional","WEIGHTS","The weights of the distributions. By default all are weighted equally.");
128 8 : keys.use("WELLTEMPERED_FACTOR");
129 8 : keys.use("SHIFT_TO_ZERO");
130 8 : keys.use("NORMALIZE");
131 8 : }
132 :
133 :
134 6 : TD_ExponentiallyModifiedGaussian::TD_ExponentiallyModifiedGaussian(const ActionOptions& ao):
135 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
136 12 : centers_(0),
137 6 : sigmas_(0),
138 6 : lambdas_(0),
139 6 : weights_(0),
140 12 : ncenters_(0) {
141 9 : for(unsigned int i=1;; i++) {
142 : std::vector<double> tmp_center;
143 30 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {
144 : break;
145 : }
146 9 : centers_.push_back(tmp_center);
147 9 : }
148 9 : for(unsigned int i=1;; i++) {
149 : std::vector<double> tmp_sigma;
150 30 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {
151 : break;
152 : }
153 20 : for(unsigned int k=0; k<tmp_sigma.size(); k++) {
154 11 : if(tmp_sigma[k]<=0.0) {
155 0 : plumed_merror(getName()+": the values given in SIGMA should be positive");
156 : }
157 : }
158 9 : sigmas_.push_back(tmp_sigma);
159 9 : }
160 9 : for(unsigned int i=1;; i++) {
161 : std::vector<double> tmp_lambda;
162 30 : if(!parseNumberedVector("LAMBDA",i,tmp_lambda) ) {
163 : break;
164 : }
165 20 : for(unsigned int k=0; k<tmp_lambda.size(); k++) {
166 11 : if(tmp_lambda[k]<=0.0) {
167 0 : plumed_merror(getName()+": the values given in LAMBDA should be positive");
168 : }
169 : }
170 9 : lambdas_.push_back(tmp_lambda);
171 9 : }
172 : //
173 6 : if(centers_.size()==0) {
174 0 : plumed_merror(getName()+": CENTER keywords seem to be missing. Note that numbered keywords start at CENTER1.");
175 : }
176 : //
177 6 : if(centers_.size()!=sigmas_.size() || centers_.size()!=lambdas_.size() ) {
178 0 : plumed_merror(getName()+": there has to be an equal amount of CENTER, SIGMA, and LAMBDA keywords");
179 : }
180 : //
181 6 : setDimension(centers_[0].size());
182 6 : ncenters_ = centers_.size();
183 : //
184 : // check centers and sigmas
185 15 : for(unsigned int i=0; i<ncenters_; i++) {
186 9 : if(centers_[i].size()!=getDimension()) {
187 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
188 : }
189 9 : if(sigmas_[i].size()!=getDimension()) {
190 0 : plumed_merror(getName()+": one of the SIGMA keyword does not match the given dimension");
191 : }
192 9 : if(lambdas_[i].size()!=getDimension()) {
193 0 : plumed_merror(getName()+": one of the LAMBDA keyword does not match the given dimension");
194 : }
195 : }
196 : //
197 12 : parseVector("WEIGHTS",weights_);
198 6 : if(weights_.size()==0) {
199 4 : weights_.assign(centers_.size(),1.0);
200 : }
201 6 : if(centers_.size()!=weights_.size()) {
202 0 : plumed_merror(getName()+": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
203 : }
204 : //
205 : double sum_weights=0.0;
206 15 : for(unsigned int i=0; i<weights_.size(); i++) {
207 9 : sum_weights+=weights_[i];
208 : }
209 15 : for(unsigned int i=0; i<weights_.size(); i++) {
210 9 : weights_[i]/=sum_weights;
211 : }
212 : //
213 6 : checkRead();
214 6 : }
215 :
216 :
217 11206 : double TD_ExponentiallyModifiedGaussian::getValue(const std::vector<double>& argument) const {
218 : double value=0.0;
219 33015 : for(unsigned int i=0; i<ncenters_; i++) {
220 21809 : value+=weights_[i]*ExponentiallyModifiedGaussianDiagonal(argument,centers_[i],sigmas_[i],lambdas_[i]);
221 : }
222 11206 : return value;
223 : }
224 :
225 :
226 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 {
227 : double value = 1.0;
228 64020 : for(unsigned int k=0; k<argument.size(); k++) {
229 42211 : double arg1 = 0.5*lambda[k]*(2.0*center[k]+lambda[k]*sigma[k]*sigma[k]-2.0*argument[k]);
230 42211 : double arg2 = (center[k]+lambda[k]*sigma[k]*sigma[k]-argument[k])/(sqrt(2.0)*sigma[k]);
231 42211 : value *= 0.5*lambda[k]*exp(arg1)*erfc(arg2);
232 : }
233 21809 : return value;
234 : }
235 :
236 :
237 :
238 : }
239 : }
|
