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