Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 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 "GridVessel.h"
23 : #include "vesselbase/ActionWithVessel.h"
24 : #include "tools/Tools.h"
25 :
26 : namespace PLMD {
27 : namespace gridtools {
28 :
29 36 : void GridVessel::registerKeywords( Keywords& keys ) {
30 36 : AveragingVessel::registerKeywords( keys );
31 36 : keys.add("compulsory","COMPONENTS","the names of the components in the vector");
32 36 : keys.add("compulsory","COORDINATES","the names of the coordinates of the grid");
33 36 : keys.add("compulsory","PBC","is the grid periodic in each direction or not");
34 36 : }
35 :
36 36 : GridVessel::GridVessel( const vesselbase::VesselOptions& da ):
37 : AveragingVessel(da),
38 : bounds_set(false),
39 : cube_units(1.0),
40 : noderiv(false),
41 36 : npoints(0)
42 : {
43 36 : std::vector<std::string> compnames; parseVector("COMPONENTS",compnames);
44 72 : std::vector<std::string> coordnames; parseVector("COORDINATES",coordnames);
45 36 : dimension=coordnames.size();
46 72 : std::vector<std::string> spbc( dimension ); parseVector("PBC",spbc);
47 36 : str_min.resize( dimension); str_max.resize( dimension ); stride.resize( dimension );
48 36 : max.resize( dimension ); dx.resize( dimension ); nbin.resize( dimension ); min.resize( dimension );
49 :
50 :
51 36 : unsigned n=0; nper=compnames.size()*( 1 + coordnames.size() );
52 36 : arg_names.resize( coordnames.size() + compnames.size()*( 1 + coordnames.size() ) );
53 36 : for(unsigned i=0; i<coordnames.size(); ++i) { arg_names[n] = coordnames[i]; n++; }
54 72 : for(unsigned i=0; i<compnames.size(); ++i) {
55 36 : arg_names[n]=compnames[i]; n++;
56 36 : for(unsigned j=0; j<coordnames.size(); ++j) { arg_names[n] = "d" + compnames[i] + "_" + coordnames[j]; n++; }
57 : }
58 :
59 36 : pbc.resize( dimension );
60 91 : for(unsigned i=0; i<dimension; ++i) {
61 55 : if( spbc[i]=="F" ) pbc[i]=false;
62 20 : else if( spbc[i]=="T" ) pbc[i]=true;
63 0 : else plumed_error();
64 36 : }
65 36 : }
66 :
67 6 : void GridVessel::setNoDerivatives() {
68 6 : nper = ( nper/(1+dimension) ); noderiv=true;
69 12 : std::vector<std::string> tnames( dimension ), cnames(nper);
70 6 : for(unsigned i=0; i<dimension; ++i) tnames[i]=arg_names[i];
71 6 : unsigned k=dimension; for(unsigned i=0; i<nper; ++i) { cnames[i]=arg_names[k]; k+=(1+dimension); }
72 6 : arg_names.resize( dimension + nper );
73 6 : for(unsigned i=0; i<dimension; ++i) arg_names[i]=tnames[i];
74 12 : for(unsigned i=0; i<nper; ++i) arg_names[dimension+i]=cnames[i];
75 6 : }
76 :
77 56 : void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
78 : const std::vector<unsigned>& binsin, const std::vector<double>& spacing ) {
79 : plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
80 56 : plumed_assert( (spacing.size()==dimension || binsin.size()==dimension) );
81 :
82 56 : npoints=1; bounds_set=true;
83 146 : for(unsigned i=0; i<dimension; ++i) {
84 90 : str_min[i]=smin[i]; str_max[i]=smax[i];
85 90 : Tools::convert( str_min[i], min[i] );
86 90 : Tools::convert( str_max[i], max[i] );
87 90 : if( spacing.size()==dimension && binsin.size()==dimension ) {
88 27 : double range = max[i] - min[i]; unsigned spc = std::floor( range / spacing[i]);
89 : // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
90 27 : if( fabs( binsin[i]*spacing[i]-range )>epsilon ) spc += 1;
91 27 : if( spc>binsin[i] ) nbin[i]=spc; else nbin[i]=binsin[i];
92 63 : } else if( binsin.size()==dimension ) nbin[i]=binsin[i];
93 0 : else if( spacing.size()==dimension ) nbin[i] = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
94 0 : else plumed_error();
95 90 : dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
96 90 : if( !pbc[i] ) { max[i] +=dx[i]; nbin[i]+=1; }
97 90 : stride[i]=npoints;
98 90 : npoints*=nbin[i];
99 : }
100 56 : resize(); // Always resize after setting new bounds as grid size may have have changed
101 56 : }
102 :
103 36 : std::string GridVessel::description() {
104 36 : if( !bounds_set ) return "";
105 :
106 48 : std::string des="grid of "; std::string num;
107 36 : for(unsigned i=0; i<dimension-1; ++i) {
108 12 : Tools::convert( nbin[i], num );
109 12 : des += num + " X ";
110 : }
111 24 : Tools::convert( nbin[dimension-1], num );
112 24 : des += num + " equally spaced points between (";
113 24 : for(unsigned i=0; i<dimension-1; ++i) des += str_min[i] + ",";
114 24 : Tools::convert( nbin[dimension-1], num );
115 24 : des += str_min[dimension-1] + ") and (";
116 24 : for(unsigned i=0; i<dimension-1; ++i) des += str_max[i] + ",";
117 24 : des += str_max[dimension-1] + ")";
118 48 : return des;
119 : }
120 :
121 109 : void GridVessel::resize() {
122 109 : plumed_massert( nper>0, "Number of datapoints at each grid point has not been set");
123 109 : resizeBuffer( getNumberOfBufferPoints()*nper ); setDataSize( npoints*nper );
124 109 : if( active.size()!=npoints) active.resize( npoints, true );
125 109 : }
126 :
127 15637495 : unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
128 : plumed_dbg_assert( bounds_set && indices.size()==dimension );
129 : // indices are flattended using a column-major order
130 15637495 : unsigned index=indices[dimension-1];
131 46883540 : for(unsigned i=dimension-1; i>0; --i) {
132 31250000 : index=index*nbin[i-1]+indices[i-1];
133 : }
134 15633540 : return index;
135 : }
136 :
137 15398 : void GridVessel::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
138 : plumed_dbg_assert( bounds_set && point.size()==dimension && indices.size()==dimension );
139 60844 : for(unsigned i=0; i<dimension; ++i) {
140 45425 : indices[i]=std::floor( (point[i] - min[i])/dx[i] );
141 45442 : if( pbc[i] ) indices[i]=indices[i]%nbin[i];
142 4881 : else if( indices[i]>nbin[i] ) plumed_merror("point is outside grid range");
143 : }
144 15419 : }
145 :
146 7516 : unsigned GridVessel::getIndex( const std::vector<double>& point ) const {
147 : plumed_dbg_assert( bounds_set && point.size()==dimension );
148 7516 : std::vector<unsigned> indices(dimension); getIndices( point, indices );
149 7516 : return getIndex( indices );
150 : }
151 :
152 76468697 : void GridVessel::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
153 76468697 : unsigned kk=index;
154 76468697 : indices[0]=index%nnbin[0];
155 152648942 : for(unsigned i=1; i<dimension-1; ++i) {
156 76231686 : kk=(kk-indices[i-1])/nnbin[i-1];
157 76510743 : indices[i]=kk%nnbin[i];
158 : }
159 76417256 : if(dimension>=2) { // I think this is wrong
160 76436052 : indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
161 : }
162 :
163 76331280 : }
164 :
165 10765265 : void GridVessel::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
166 10765265 : convertIndexToIndices( index, nbin, indices );
167 10764566 : }
168 :
169 10664906 : void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
170 : plumed_dbg_assert( x.size()==dimension && ipoint<npoints );
171 10664906 : std::vector<unsigned> tindices( dimension ); getIndices( ipoint, tindices );
172 10667313 : for(unsigned i=0; i<dimension; ++i) x[i] = min[i] + dx[i]*tindices[i];
173 10665087 : }
174 :
175 7510 : void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const {
176 7510 : mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
177 :
178 7511 : std::vector<unsigned> tmp_indices( dimension );
179 15022 : std::vector<unsigned> my_indices( dimension );
180 7511 : getIndices( mybox, my_indices );
181 66389 : for(unsigned i=0; i<mysneigh.size(); ++i) {
182 58868 : unsigned tmp=i;
183 234657 : for(unsigned j=0; j<dimension; ++j) {
184 175774 : unsigned i0=tmp%2+my_indices[j]; tmp/=2;
185 175807 : if(!pbc[j] && i0==nbin[j]) getAction()->error("Extrapolating function on grid");
186 175804 : if( pbc[j] && i0==nbin[j]) i0=0;
187 175791 : tmp_indices[j]=i0;
188 : }
189 58883 : mysneigh[i]=getIndex( tmp_indices );
190 58880 : plumed_massert( active[mysneigh[i]], "inactive grid point required for splines");
191 7511 : }
192 7511 : }
193 :
194 46612 : double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
195 46612 : plumed_assert( bounds_set && ipoint<npoints && jelement<nper && active[ipoint] );
196 46583 : return getDataElement( nper*ipoint + jelement );
197 : }
198 :
199 0 : void GridVessel::setGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ) {
200 : plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
201 0 : setDataElement( nper*ipoint + jelement, value );
202 0 : }
203 :
204 5 : void GridVessel::addToGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ) {
205 : plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
206 5 : addDataElement( nper*ipoint + jelement, value );
207 5 : }
208 :
209 9204 : void GridVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
210 : plumed_dbg_assert( myvals.getNumberOfValues()==(nper+1) );
211 9204 : for(unsigned i=0; i<nper; ++i) buffer[bufstart + nper*current + i] += myvals.get(i+1);
212 9279 : }
213 :
214 0 : double GridVessel::getGridElement( const std::vector<unsigned>& indices, const unsigned& jelement ) const {
215 0 : return getGridElement( getIndex( indices ), jelement );
216 : }
217 :
218 0 : void GridVessel::setGridElement( const std::vector<unsigned>& indices, const unsigned& jelement, const double& value ) {
219 0 : setGridElement( getIndex( indices ), jelement, value );
220 0 : }
221 :
222 4837 : std::vector<std::string> GridVessel::getMin() const {
223 4837 : return str_min;
224 : }
225 :
226 4866 : std::vector<std::string> GridVessel::getMax() const {
227 4866 : return str_max;
228 : }
229 :
230 5259 : std::vector<unsigned> GridVessel::getNbin() const {
231 : plumed_dbg_assert( bounds_set );
232 5259 : std::vector<unsigned> ngrid( dimension );
233 12182 : for(unsigned i=0; i<dimension; ++i) {
234 6921 : if( !pbc[i] ) ngrid[i]=nbin[i] - 1;
235 3492 : else ngrid[i]=nbin[i];
236 : }
237 5261 : return ngrid;
238 : }
239 :
240 15996 : void GridVessel::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
241 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
242 : plumed_dbg_assert( bounds_set && nneigh.size()==dimension );
243 :
244 15996 : std::vector<unsigned> indices( dimension );
245 16000 : for(unsigned i=0; i<dimension; ++i) indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
246 15997 : getNeighbors( indices, nneigh, num_neighbors, neighbors );
247 15998 : }
248 :
249 34083 : void GridVessel::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
250 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
251 : plumed_dbg_assert( bounds_set && nneigh.size()==dimension );
252 :
253 34083 : unsigned num_neigh=1; std::vector<unsigned> small_bin( dimension );
254 136256 : for(unsigned i=0; i<dimension; ++i) {
255 102170 : small_bin[i]=(2*nneigh[i]+1);
256 102170 : num_neigh *=small_bin[i];
257 : }
258 34086 : if( neighbors.size()!=num_neigh ) neighbors.resize( num_neigh );
259 :
260 34086 : num_neighbors=0;
261 68172 : std::vector<unsigned> s_indices(dimension), t_indices(dimension);
262 65813929 : for(unsigned index=0; index<num_neigh; ++index) {
263 65779843 : bool found=true;
264 65779843 : convertIndexToIndices( index, small_bin, s_indices );
265 262809041 : for(unsigned i=0; i<dimension; ++i) {
266 197053097 : int i0=s_indices[i]-nneigh[i]+indices[i];
267 197097486 : if(!pbc[i] && i0<0) found=false;
268 197368259 : if(!pbc[i] && i0>=nbin[i]) found=false;
269 197149502 : if( pbc[i] && i0<0) i0=nbin[i]-(-i0)%nbin[i];
270 197142968 : if( pbc[i] && i0>=nbin[i]) i0%=nbin[i];
271 197044415 : t_indices[i]=static_cast<unsigned>(i0);
272 : }
273 65755944 : if( found ) {
274 15491747 : neighbors[num_neighbors]=getIndex( t_indices );
275 15490356 : num_neighbors++;
276 : }
277 34086 : }
278 34086 : }
279 :
280 11 : void GridVessel::setCubeUnits( const double& units ) {
281 11 : cube_units=units;
282 11 : }
283 :
284 8 : double GridVessel::getCubeUnits() const {
285 8 : return cube_units;
286 : }
287 :
288 7 : std::string GridVessel::getInputString() const {
289 7 : std::string mstring="COORDINATES="+arg_names[0];
290 7 : for(unsigned i=1; i<dimension; ++i) mstring+="," + arg_names[i];
291 7 : mstring += " PBC=";
292 7 : if( pbc[0] ) mstring +="T";
293 5 : else mstring +="F";
294 11 : for(unsigned i=1; i<dimension; ++i) {
295 4 : if( pbc[i] ) mstring +=",T";
296 4 : else mstring +=",F";
297 : }
298 7 : return mstring;
299 : }
300 :
301 7510 : double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const {
302 : plumed_dbg_assert( der.size()==dimension && !noderiv && ind<getNumberOfComponents() );
303 :
304 7510 : double X,X2,X3,value=0; der.assign(der.size(),0.0);
305 7511 : std::vector<double> fd(dimension);
306 15022 : std::vector<double> C(dimension);
307 15022 : std::vector<double> D(dimension);
308 15021 : std::vector<double> dder(dimension);
309 :
310 15022 : std::vector<unsigned> nindices(dimension);
311 15022 : std::vector<unsigned> indices(dimension); getIndices( x, indices );
312 15022 : std::vector<unsigned> neigh; getSplineNeighbors( getIndex(indices), neigh );
313 15021 : std::vector<double> xfloor(dimension); getGridPointCoordinates( getIndex(x), xfloor );
314 :
315 : // loop over neighbors
316 66390 : for(unsigned int ipoint=0; ipoint<neigh.size(); ++ipoint) {
317 58871 : double grid=getGridElement(neigh[ipoint], ind*(1+dimension) );
318 58913 : for(unsigned j=0; j<dimension; ++j) dder[j] = getGridElement( neigh[ipoint], ind*(1+dimension) + 1 + j );
319 :
320 58879 : getIndices( neigh[ipoint], nindices );
321 58823 : double ff=1.0;
322 :
323 234607 : for(unsigned j=0; j<dimension; ++j) {
324 175742 : int x0=1;
325 175742 : if(nindices[j]==indices[j]) x0=0;
326 175801 : double ddx=dx[j];
327 175787 : X=fabs((x[j]-xfloor[j])/ddx-(double)x0);
328 175613 : X2=X*X;
329 175613 : X3=X2*X;
330 : double yy;
331 175613 : if(fabs(grid)<0.0000001) yy=0.0;
332 161521 : else yy=-dder[j]/grid;
333 175823 : C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*ddx;
334 175857 : D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*ddx;
335 175837 : D[j]*=(x0?-1.0:1.0)/ddx;
336 175849 : ff*=C[j];
337 : }
338 234650 : for(unsigned j=0; j<dimension; ++j) {
339 175794 : fd[j]=D[j];
340 175901 : for(unsigned i=0; i<dimension; ++i) if(i!=j) fd[j]*=C[i];
341 : }
342 58856 : value+=grid*ff;
343 58856 : for(unsigned j=0; j<dimension; ++j) der[j]+=grid*fd[j];
344 : }
345 15021 : return value;
346 : }
347 :
348 2 : void GridVessel::activateThesePoints( const std::vector<bool>& to_activate ) {
349 : plumed_dbg_assert( to_activate.size()==npoints );
350 2 : for(unsigned i=0; i<npoints; ++i) active[i]=to_activate[i];
351 2 : }
352 :
353 : }
354 2523 : }
355 :
|
