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 : #include "GridIntegrationWeights.h"
25 : #include "VesTools.h"
26 : #include "GridLinearInterpolation.h"
27 :
28 : #include "core/ActionRegister.h"
29 : #include "tools/Grid.h"
30 : #include "core/Value.h"
31 : #include "tools/File.h"
32 :
33 :
34 : namespace PLMD {
35 : namespace ves {
36 :
37 :
38 : //+PLUMEDOC VES_TARGETDIST TD_GRID
39 : /*
40 : Target distribution from an external grid file (static).
41 :
42 : Using this keyword you can use a target distribution that is read from an
43 : external grid file that is in the proper PLUMED file format. You do not to
44 : give any information about the external grid file as all relevant information
45 : should be automatically detected. It is assumed that the distribution read in
46 : from the grid is a proper probability distribution, i.e. always non-negative
47 : and can be normalized.
48 :
49 : By default the target distribution from the external grid is always normalized
50 : inside the code. You can disable this normalization by using DO_NOT_NORMALIZE
51 : keyword. However, be warned that this will generally lead to the wrong
52 : behavior if the distribution from the external grid is not properly
53 : normalized to 1.
54 :
55 : If the distribution from the external grid file has for some reason
56 : negative values can you use the SHIFT keyword to shift the distribution
57 : by a given value. Another option is to use
58 : the SHIFT_TO_ZERO keyword to shift the minimum of the distribution to zero.
59 :
60 : Note that the number of grid bins used in the external grid file do not have
61 : to be the same as used in the bias or action where the target distribution is
62 : employed as the code will employ a linear (or bilinear for two dimensions)
63 : interpolation to calculate values. Currently only one or two dimensional grids
64 : are supported.
65 :
66 : It can happen that the intervals on which the target distribution is defined is
67 : larger than the intervals covered by the external grid file. In this case the
68 : default option is to consider the target distribution as continuous such that
69 : values outside the boundary of the external grid file are the same as at
70 : the boundary. This can be changed by using the ZERO_OUTSIDE keyword which
71 : will make values outside to be taken as zero.
72 :
73 : \par Examples
74 :
75 : Generally you only need to provide the the filename of the external grid
76 : file.
77 : \plumedfile
78 : td: TD_GRID FILE=input-grid.data
79 : \endplumedfile
80 :
81 : */
82 : //+ENDPLUMEDOC
83 :
84 :
85 : class TD_Grid : public TargetDistribution {
86 : Grid* distGrid_;
87 : std::vector<double> minima_;
88 : std::vector<double> maxima_;
89 : std::vector<bool> periodic_;
90 : bool zero_outside_;
91 : double shift_;
92 : public:
93 : static void registerKeywords( Keywords&);
94 : explicit TD_Grid(const ActionOptions& ao);
95 : ~TD_Grid();
96 : double getValue(const std::vector<double>&) const ;
97 : };
98 :
99 :
100 4830 : PLUMED_REGISTER_ACTION(TD_Grid,"TD_GRID")
101 :
102 :
103 10 : void TD_Grid::registerKeywords(Keywords& keys) {
104 10 : TargetDistribution::registerKeywords(keys);
105 10 : keys.add("compulsory","FILE","The name of the external grid file to be used as a target distribution.");
106 10 : keys.add("optional","SHIFT","Shift the grid read in by some constant value. Due to normalization the final shift in the target distribution will generally not be the same as the value given here");
107 10 : keys.addFlag("ZERO_OUTSIDE",false,"By default the target distribution is continuous such that values outside the boundary of the external grid file are the same as at the boundary. This can be changed by using this flag which will make values outside to be taken as zero.");
108 10 : keys.addFlag("DO_NOT_NORMALIZE",false,"By default the target distribution from the external grid is always normalized inside the code. You can use this flag to disable this normalization. However, be warned that this will generally lead to the wrong behavior if the distribution from the external grid is not properly normalized to 1.");
109 10 : keys.use("WELLTEMPERED_FACTOR");
110 10 : keys.use("SHIFT_TO_ZERO");
111 10 : }
112 :
113 27 : TD_Grid::~TD_Grid() {
114 9 : if(distGrid_!=NULL) {
115 9 : delete distGrid_;
116 : }
117 18 : }
118 :
119 :
120 9 : TD_Grid::TD_Grid(const ActionOptions& ao):
121 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
122 : distGrid_(NULL),
123 : minima_(0),
124 : maxima_(0),
125 : zero_outside_(false),
126 9 : shift_(0.0)
127 : {
128 9 : std::string filename;
129 9 : parse("FILE",filename);
130 9 : parse("SHIFT",shift_);
131 9 : if(shift_!=0.0) {
132 1 : if(isTargetDistGridShiftedToZero()) {plumed_merror(getName() + ": using both SHIFT and SHIFT_TO_ZERO is not allowed.");}
133 : }
134 9 : parseFlag("ZERO_OUTSIDE",zero_outside_);
135 :
136 9 : bool do_not_normalize=false;
137 9 : parseFlag("DO_NOT_NORMALIZE",do_not_normalize);
138 9 : if(do_not_normalize && shift_!=0.0) {plumed_merror(getName() + ": using both SHIFT and DO_NOT_NORMALIZE is not allowed.");}
139 9 : if(!do_not_normalize) {setForcedNormalization();}
140 :
141 9 : checkRead();
142 :
143 18 : std::string gridlabel;
144 18 : std::vector<std::string> arglabels;
145 18 : std::vector<std::string> argmin;
146 18 : std::vector<std::string> argmax;
147 18 : std::vector<bool> argperiodic;
148 18 : std::vector<unsigned int> argnbins;
149 9 : bool has_deriv = false;
150 9 : unsigned int nargs = VesTools::getGridFileInfo(filename,gridlabel,arglabels,argmin,argmax,argperiodic,argnbins,has_deriv);
151 9 : if(nargs==0) {
152 0 : plumed_merror(getName() + ": problem in parsing information from grid file");
153 : }
154 :
155 9 : setDimension(arglabels.size());
156 18 : std::vector<Value*> arguments(arglabels.size());
157 22 : for(unsigned int i=0; i < arglabels.size(); i++) {
158 13 : arguments[i]= new Value(NULL,arglabels[i],false);
159 13 : if(argperiodic[i]) {
160 0 : arguments[i]->setDomain(argmin[i],argmax[i]);
161 : }
162 : else {
163 13 : arguments[i]->setNotPeriodic();
164 : }
165 : }
166 :
167 18 : IFile gridfile; gridfile.open(filename);
168 9 : if(has_deriv) {
169 0 : distGrid_=Grid::create(gridlabel,arguments,gridfile,false,true,true);
170 : }
171 : else {
172 9 : distGrid_=Grid::create(gridlabel,arguments,gridfile,false,false,false);
173 : }
174 9 : gridfile.close();
175 :
176 :
177 9 : minima_.resize(getDimension());
178 9 : maxima_.resize(getDimension());
179 9 : periodic_.resize(getDimension());
180 22 : for (unsigned int i=0; i < getDimension(); i++) {
181 13 : Tools::convert(distGrid_->getMin()[i],minima_[i]);
182 13 : Tools::convert(distGrid_->getMax()[i],maxima_[i]);
183 13 : periodic_[i] = argperiodic[i];
184 13 : if(periodic_[i]) {maxima_[i]-=distGrid_->getDx()[i];}
185 : }
186 :
187 9 : for(unsigned int i=0; i < arguments.size(); i++) {delete arguments[i];}
188 18 : arguments.clear();
189 9 : }
190 :
191 :
192 49369 : double TD_Grid::getValue(const std::vector<double>& argument) const {
193 49369 : double outside = 0.0;
194 49369 : std::vector<double> arg = argument;
195 145566 : for(unsigned int k=0; k<getDimension(); k++) {
196 96533 : if(zero_outside_ && (argument[k] < minima_[k] || argument[k] > maxima_[k])) {
197 336 : return outside;
198 : }
199 96197 : else if(argument[k] < minima_[k]) {
200 168 : arg[k] = minima_[k];
201 : }
202 96029 : else if(argument[k] > maxima_[k]) {
203 168 : arg[k] =maxima_[k];
204 : }
205 : }
206 49033 : return GridLinearInterpolation::getGridValueWithLinearInterpolation(distGrid_,arg)+shift_;
207 : }
208 :
209 :
210 : }
211 4821 : }
|
