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 6461 : PLUMED_REGISTER_ACTION(TD_Grid,"TD_GRID")
101 :
102 :
103 10 : void TD_Grid::registerKeywords(Keywords& keys) {
104 10 : TargetDistribution::registerKeywords(keys);
105 40 : keys.add("compulsory","FILE","The name of the external grid file to be used as a target distribution.");
106 40 : 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 30 : 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 30 : 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 20 : keys.use("WELLTEMPERED_FACTOR");
110 20 : 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 18 : shift_(0.0)
127 : {
128 : std::string filename;
129 18 : parse("FILE",filename);
130 18 : 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 18 : parseFlag("ZERO_OUTSIDE",zero_outside_);
135 :
136 9 : bool do_not_normalize=false;
137 18 : 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 : std::string gridlabel;
144 9 : std::vector<std::string> arglabels;
145 9 : std::vector<std::string> argmin;
146 9 : std::vector<std::string> argmax;
147 : std::vector<bool> argperiodic;
148 : 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 9 : std::vector<Value*> arguments(arglabels.size());
157 57 : for(unsigned int i=0; i < arglabels.size(); i++) {
158 26 : 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 35 : for (unsigned int i=0; i < getDimension(); i++) {
181 39 : Tools::convert(distGrid_->getMin()[i],minima_[i]);
182 26 : Tools::convert(distGrid_->getMax()[i],maxima_[i]);
183 : periodic_[i] = argperiodic[i];
184 13 : if(periodic_[i]) {maxima_[i]-=distGrid_->getDx()[i];}
185 : }
186 :
187 57 : for(unsigned int i=0; i < arguments.size(); i++) {delete arguments[i];}
188 : arguments.clear();
189 9 : }
190 :
191 :
192 49369 : double TD_Grid::getValue(const std::vector<double>& argument) const {
193 : double outside = 0.0;
194 49369 : std::vector<double> arg = argument;
195 241763 : for(unsigned int k=0; k<getDimension(); k++) {
196 308429 : if(zero_outside_ && (argument[k] < minima_[k] || argument[k] > maxima_[k])) {
197 : return outside;
198 : }
199 288591 : 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 4839 : }
|
