Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed 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 : plumed 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 plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "BiasRepresentation.h"
23 : #include "core/Value.h"
24 : #include "Communicator.h"
25 : #include <iostream>
26 : #include "KernelFunctions.h"
27 : #include "File.h"
28 : #include "Grid.h"
29 :
30 : namespace PLMD {
31 :
32 : /// the constructor here
33 4 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc ):hasgrid(false),rescaledToBias(false),mycomm(cc) {
34 4 : lowI_=0.0;
35 4 : uppI_=0.0;
36 4 : doInt_=false;
37 4 : ndim=tmpvalues.size();
38 12 : for(int i=0; i<ndim; i++) {
39 8 : values.push_back(tmpvalues[i]);
40 8 : names.push_back(values[i]->getName());
41 : }
42 4 : }
43 :
44 : /// overload the constructor: add the sigma at constructor time
45 1 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc, const std::vector<double> & sigma ):
46 1 : hasgrid(false), rescaledToBias(false), histosigma(sigma),mycomm(cc) {
47 1 : lowI_=0.0;
48 1 : uppI_=0.0;
49 1 : doInt_=false;
50 1 : ndim=tmpvalues.size();
51 3 : for(int i=0; i<ndim; i++) {
52 2 : values.push_back(tmpvalues[i]);
53 2 : names.push_back(values[i]->getName());
54 : }
55 1 : }
56 :
57 : /// overload the constructor: add the grid at constructor time
58 8 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc, const std::vector<std::string> & gmin, const std::vector<std::string> & gmax,
59 8 : const std::vector<unsigned> & nbin, bool doInt, double lowI, double uppI):
60 8 : hasgrid(false), rescaledToBias(false), mycomm(cc) {
61 8 : ndim=tmpvalues.size();
62 23 : for(int i=0; i<ndim; i++) {
63 15 : values.push_back(tmpvalues[i]);
64 15 : names.push_back(values[i]->getName());
65 : }
66 8 : doInt_=doInt;
67 8 : lowI_=lowI;
68 8 : uppI_=uppI;
69 : // initialize the grid
70 8 : addGrid(gmin,gmax,nbin);
71 8 : }
72 :
73 : /// overload the constructor with some external sigmas: needed for histogram
74 1 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc, const std::vector<std::string> & gmin, const std::vector<std::string> & gmax,
75 1 : const std::vector<unsigned> & nbin, const std::vector<double> & sigma):
76 1 : hasgrid(false), rescaledToBias(false),histosigma(sigma),mycomm(cc) {
77 1 : lowI_=0.0;
78 1 : uppI_=0.0;
79 1 : doInt_=false;
80 1 : ndim=tmpvalues.size();
81 3 : for(int i=0; i<ndim; i++) {
82 2 : values.push_back(tmpvalues[i]);
83 2 : names.push_back(values[i]->getName());
84 : }
85 : // initialize the grid
86 1 : addGrid(gmin,gmax,nbin);
87 1 : }
88 :
89 9 : void BiasRepresentation::addGrid(const std::vector<std::string> & gmin, const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin ) {
90 9 : plumed_massert(hills.size()==0,"you can set the grid before loading the hills");
91 9 : plumed_massert(hasgrid==false,"to build the grid you should not having the grid in this bias representation");
92 : std::string ss;
93 : ss="file.free";
94 : std::vector<Value*> vv;
95 26 : for(unsigned i=0; i<values.size(); i++) {
96 17 : vv.push_back(values[i]);
97 : }
98 18 : BiasGrid_=Tools::make_unique<Grid>(ss,vv,gmin,gmax,nbin,false,true);
99 9 : hasgrid=true;
100 9 : }
101 :
102 5563 : bool BiasRepresentation::hasSigmaInInput() {
103 5563 : if(histosigma.size()==0) {
104 : return false;
105 : } else {
106 1092 : return true;
107 : }
108 : }
109 :
110 1 : void BiasRepresentation::setRescaledToBias(bool rescaled) {
111 1 : plumed_massert(hills.size()==0,"you can set the rescaling function only before loading hills");
112 1 : rescaledToBias=rescaled;
113 1 : }
114 :
115 0 : const bool & BiasRepresentation::isRescaledToBias() {
116 0 : return rescaledToBias;
117 : }
118 :
119 0 : unsigned BiasRepresentation::getNumberOfDimensions() {
120 0 : return values.size();
121 : }
122 :
123 0 : std::vector<std::string> BiasRepresentation::getNames() {
124 0 : return names;
125 : }
126 :
127 0 : const std::string & BiasRepresentation::getName(unsigned i) {
128 0 : return names[i];
129 : }
130 :
131 0 : const std::vector<Value*>& BiasRepresentation::getPtrToValues() {
132 0 : return values;
133 : }
134 :
135 0 : Value* BiasRepresentation::getPtrToValue(unsigned i) {
136 0 : return values[i];
137 : }
138 :
139 1092 : std::unique_ptr<KernelFunctions> BiasRepresentation::readFromPoint(IFile *ifile) {
140 1092 : std::vector<double> cc( names.size() );
141 3276 : for(unsigned i=0; i<names.size(); ++i) {
142 2184 : ifile->scanField(names[i],cc[i]);
143 : }
144 1092 : double h=1.0;
145 2184 : return Tools::make_unique<KernelFunctions>(cc,histosigma,"stretched-gaussian","DIAGONAL",h);
146 : }
147 :
148 5563 : void BiasRepresentation::pushKernel( IFile *ifile ) {
149 5563 : std::unique_ptr<KernelFunctions> kk;
150 : // here below the reading of the kernel is completely hidden
151 5563 : if(histosigma.size()==0) {
152 4471 : ifile->allowIgnoredFields();
153 8942 : kk=KernelFunctions::read(ifile,true,names);
154 : } else {
155 : // when doing histogram assume gaussian with a given diagonal sigma
156 : // and neglect all the rest
157 2184 : kk=readFromPoint(ifile);
158 : }
159 : // the bias factor is not something about the kernels but
160 : // must be stored to keep the bias/free energy duality
161 : std::string dummy;
162 : double dummyd;
163 11126 : if(ifile->FieldExist("biasf")) {
164 5563 : ifile->scanField("biasf",dummy);
165 5563 : Tools::convert(dummy,dummyd);
166 : } else {
167 0 : dummyd=1.0;
168 : }
169 5563 : biasf.push_back(dummyd);
170 : // the domain does not pertain to the kernel but to the values here defined
171 : std::string mins,maxs,minv,maxv,mini,maxi;
172 : mins="min_";
173 : maxs="max_";
174 16688 : for(int i=0 ; i<ndim; i++) {
175 11125 : if(values[i]->isPeriodic()) {
176 22216 : ifile->scanField(mins+names[i],minv);
177 22216 : ifile->scanField(maxs+names[i],maxv);
178 : // verify that the domain is correct
179 11108 : values[i]->getDomain(mini,maxi);
180 11108 : plumed_massert(mini==minv,"the input periodicity in hills and in value definition does not match" );
181 11108 : plumed_massert(maxi==maxv,"the input periodicity in hills and in value definition does not match" );
182 : }
183 : }
184 : // if grid is defined then it should be added on the grid
185 : //cerr<<"now with "<<hills.size()<<endl;
186 5563 : if(hasgrid) {
187 : std::vector<unsigned> nneighb;
188 3375 : if(doInt_&&(kk->getCenter()[0]+kk->getContinuousSupport()[0] > uppI_ || kk->getCenter()[0]-kk->getContinuousSupport()[0] < lowI_ )) {
189 0 : nneighb=BiasGrid_->getNbin();
190 : } else {
191 6750 : nneighb=kk->getSupport(BiasGrid_->getDx());
192 : }
193 3375 : std::vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(kk->getCenter(),nneighb);
194 3375 : std::vector<double> der(ndim);
195 3375 : std::vector<double> xx(ndim);
196 3375 : if(mycomm.Get_size()==1) {
197 1027680 : for(unsigned i=0; i<neighbors.size(); ++i) {
198 1024305 : Grid::index_t ineigh=neighbors[i];
199 3072823 : for(int j=0; j<ndim; ++j) {
200 2048518 : der[j]=0.0;
201 : }
202 1024305 : BiasGrid_->getPoint(ineigh,xx);
203 : // assign xx to a new vector of values
204 3072823 : for(int j=0; j<ndim; ++j) {
205 2048518 : values[j]->set(xx[j]);
206 : }
207 : double bias;
208 1024305 : if(doInt_) {
209 0 : bias=kk->evaluate(values,der,true,doInt_,lowI_,uppI_);
210 : } else {
211 1024305 : bias=kk->evaluate(values,der,true);
212 : }
213 1024305 : if(rescaledToBias) {
214 45234 : double f=(biasf.back()-1.)/(biasf.back());
215 45234 : bias*=f;
216 135702 : for(int j=0; j<ndim; ++j) {
217 90468 : der[j]*=f;
218 : }
219 : }
220 1024305 : BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
221 : }
222 : } else {
223 0 : unsigned stride=mycomm.Get_size();
224 0 : unsigned rank=mycomm.Get_rank();
225 0 : std::vector<double> allder(ndim*neighbors.size(),0.0);
226 0 : std::vector<double> allbias(neighbors.size(),0.0);
227 0 : std::vector<double> tmpder(ndim);
228 0 : for(unsigned i=rank; i<neighbors.size(); i+=stride) {
229 0 : Grid::index_t ineigh=neighbors[i];
230 0 : BiasGrid_->getPoint(ineigh,xx);
231 0 : for(int j=0; j<ndim; ++j) {
232 0 : values[j]->set(xx[j]);
233 : }
234 0 : if(doInt_) {
235 0 : allbias[i]=kk->evaluate(values,der,true,doInt_,lowI_,uppI_);
236 : } else {
237 0 : allbias[i]=kk->evaluate(values,der,true);
238 : }
239 0 : if(rescaledToBias) {
240 0 : double f=(biasf.back()-1.)/(biasf.back());
241 0 : allbias[i]*=f;
242 0 : for(int j=0; j<ndim; ++j) {
243 0 : tmpder[j]*=f;
244 : }
245 : }
246 : // this solution with the temporary vector is rather bad, probably better to take
247 : // a pointer of double as it was in old gaussian
248 0 : for(int j=0; j<ndim; ++j) {
249 0 : allder[ndim*i+j]=tmpder[j];
250 0 : tmpder[j]=0.;
251 : }
252 : }
253 0 : mycomm.Sum(allbias);
254 0 : mycomm.Sum(allder);
255 0 : for(unsigned i=0; i<neighbors.size(); ++i) {
256 0 : Grid::index_t ineigh=neighbors[i];
257 0 : for(int j=0; j<ndim; ++j) {
258 0 : der[j]=allder[ndim*i+j];
259 : }
260 0 : BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
261 : }
262 : }
263 : }
264 5563 : hills.emplace_back(std::move(kk));
265 5563 : }
266 :
267 5578 : int BiasRepresentation::getNumberOfKernels() {
268 5578 : return hills.size();
269 : }
270 :
271 10 : Grid* BiasRepresentation::getGridPtr() {
272 10 : plumed_massert(hasgrid,"if you want the grid pointer then you should have defined a grid before");
273 10 : return BiasGrid_.get();
274 : }
275 :
276 5 : void BiasRepresentation::getMinMaxBin(std::vector<double> &vmin, std::vector<double> &vmax, std::vector<unsigned> &vbin) {
277 : std::vector<double> ss,cc,binsize;
278 : vmin.clear();
279 5 : vmin.resize(ndim,10.e20);
280 : vmax.clear();
281 5 : vmax.resize(ndim,-10.e20);
282 : vbin.clear();
283 5 : vbin.resize(ndim);
284 : binsize.clear();
285 5 : binsize.resize(ndim,10.e20);
286 : int ndiv=10; // adjustable parameter: division per support
287 2193 : for(unsigned i=0; i<hills.size(); i++) {
288 2188 : if(histosigma.size()!=0) {
289 546 : ss=histosigma;
290 : } else {
291 3284 : ss=hills[i]->getContinuousSupport();
292 : }
293 2188 : cc=hills[i]->getCenter();
294 6564 : for(int j=0; j<ndim; j++) {
295 4376 : double dmin=cc[j]-ss[j];
296 4376 : double dmax=cc[j]+ss[j];
297 4376 : double ddiv=ss[j]/double(ndiv);
298 4376 : if(dmin<vmin[j]) {
299 300 : vmin[j]=dmin;
300 : }
301 4376 : if(dmax>vmax[j]) {
302 236 : vmax[j]=dmax;
303 : }
304 4376 : if(ddiv<binsize[j]) {
305 34 : binsize[j]=ddiv;
306 : }
307 : }
308 : }
309 15 : for(int j=0; j<ndim; j++) {
310 : // reset to periodicity
311 10 : if(values[j]->isPeriodic()) {
312 : double minv,maxv;
313 8 : values[j]->getDomain(minv,maxv);
314 8 : if(minv>vmin[j]) {
315 0 : vmin[j]=minv;
316 : }
317 8 : if(maxv<vmax[j]) {
318 0 : vmax[j]=maxv;
319 : }
320 : }
321 10 : vbin[j]=static_cast<unsigned>(std::ceil((vmax[j]-vmin[j])/binsize[j]) );
322 : }
323 5 : }
324 :
325 0 : void BiasRepresentation::clear() {
326 0 : hills.clear();
327 : // clear the grid
328 0 : if(hasgrid) {
329 0 : BiasGrid_->clear();
330 : }
331 0 : }
332 :
333 :
334 : }
|
