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