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