Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "GridVessel.h"
23 : #include "vesselbase/ActionWithVessel.h"
24 : #include "tools/Random.h"
25 : #include "tools/Tools.h"
26 :
27 : namespace PLMD {
28 : namespace gridtools {
29 :
30 66 : void GridVessel::registerKeywords( Keywords& keys ) {
31 66 : AveragingVessel::registerKeywords( keys );
32 132 : keys.add("compulsory","TYPE","flat","how the grid points are being generated");
33 132 : keys.add("compulsory","COMPONENTS","the names of the components in the vector");
34 132 : keys.add("compulsory","COORDINATES","the names of the coordinates of the grid");
35 132 : keys.add("compulsory","PBC","is the grid periodic in each direction or not");
36 66 : }
37 :
38 67 : GridVessel::GridVessel( const vesselbase::VesselOptions& da ):
39 : AveragingVessel(da),
40 67 : bounds_set(false),
41 67 : npoints(0),
42 67 : cube_units(1.0),
43 67 : wasforced(false),
44 67 : noderiv(false) {
45 : std::string geom;
46 134 : parse("TYPE",geom);
47 67 : if( geom=="flat" ) {
48 64 : gtype=flat;
49 3 : } else if( geom=="fibonacci" ) {
50 3 : gtype=fibonacci;
51 : } else {
52 0 : plumed_merror( geom + " is invalid geometry type");
53 : }
54 : std::vector<std::string> compnames;
55 134 : parseVector("COMPONENTS",compnames);
56 : std::vector<std::string> coordnames;
57 67 : parseVector("COORDINATES",coordnames);
58 67 : if( gtype==flat ) {
59 64 : dimension=coordnames.size();
60 64 : str_min.resize( dimension);
61 64 : str_max.resize( dimension );
62 64 : stride.resize( dimension );
63 64 : max.resize( dimension );
64 64 : dx.resize( dimension );
65 64 : nbin.resize( dimension );
66 64 : min.resize( dimension );
67 3 : } else if( gtype==fibonacci ) {
68 3 : if( coordnames.size()!=3 ) {
69 0 : error("cannot generate fibonacci grid points on surface of sphere if not 3 input coordinates");
70 : }
71 3 : dimension=3;
72 : }
73 :
74 : unsigned n=0;
75 67 : nper=compnames.size()*( 1 + coordnames.size() );
76 67 : arg_names.resize( coordnames.size() + compnames.size()*( 1 + coordnames.size() ) );
77 169 : for(unsigned i=0; i<coordnames.size(); ++i) {
78 102 : arg_names[n] = coordnames[i];
79 102 : n++;
80 : }
81 135 : for(unsigned i=0; i<compnames.size(); ++i) {
82 68 : arg_names[n]=compnames[i];
83 68 : n++;
84 172 : for(unsigned j=0; j<coordnames.size(); ++j) {
85 208 : arg_names[n] = "d" + compnames[i] + "_" + coordnames[j];
86 104 : n++;
87 : }
88 : }
89 67 : pbc.resize( dimension );
90 67 : std::vector<std::string> spbc( dimension );
91 67 : parseVector("PBC",spbc);
92 169 : for(unsigned i=0; i<dimension; ++i) {
93 102 : if( spbc[i]=="F" ) {
94 : pbc[i]=false;
95 31 : } else if( spbc[i]=="T" ) {
96 : pbc[i]=true;
97 : } else {
98 0 : plumed_error();
99 : }
100 : }
101 134 : }
102 :
103 19 : void GridVessel::setNoDerivatives() {
104 19 : nper = ( nper/(1+dimension) );
105 19 : noderiv=true;
106 19 : std::vector<std::string> tnames( dimension ), cnames(nper);
107 41 : for(unsigned i=0; i<dimension; ++i) {
108 22 : tnames[i]=arg_names[i];
109 : }
110 : unsigned k=dimension;
111 38 : for(unsigned i=0; i<nper; ++i) {
112 19 : cnames[i]=arg_names[k];
113 19 : k+=(1+dimension);
114 : }
115 19 : arg_names.resize( dimension + nper );
116 41 : for(unsigned i=0; i<dimension; ++i) {
117 22 : arg_names[i]=tnames[i];
118 : }
119 38 : for(unsigned i=0; i<nper; ++i) {
120 19 : arg_names[dimension+i]=cnames[i];
121 : }
122 19 : }
123 :
124 98 : void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
125 : const std::vector<unsigned>& binsin, const std::vector<double>& spacing ) {
126 : plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
127 98 : plumed_assert( gtype==flat && (spacing.size()==dimension || binsin.size()==dimension) );
128 :
129 98 : npoints=1;
130 98 : bounds_set=true;
131 245 : for(unsigned i=0; i<dimension; ++i) {
132 147 : str_min[i]=smin[i];
133 : str_max[i]=smax[i];
134 : // GB: I ignore the result here not to break test multicolvar/rt-dens-average
135 : // where these functions were called with str_min[i] and str_max[i] as empty string
136 : // To be checked.
137 147 : Tools::convertNoexcept( str_min[i], min[i] );
138 147 : Tools::convertNoexcept( str_max[i], max[i] );
139 :
140 147 : if( spacing.size()==dimension && binsin.size()==dimension ) {
141 39 : if( spacing[i]==0 ) {
142 2 : nbin[i] = binsin[i];
143 : } else {
144 37 : double range = max[i] - min[i];
145 37 : nbin[i] = std::round( range / spacing[i]);
146 : // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
147 37 : if( nbin[i]!=binsin[i] ) {
148 0 : plumed_merror("mismatch between input spacing and input number of bins");
149 : }
150 : }
151 108 : } else if( binsin.size()==dimension ) {
152 108 : nbin[i]=binsin[i];
153 0 : } else if( spacing.size()==dimension ) {
154 0 : nbin[i] = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
155 : } else {
156 0 : plumed_error();
157 : }
158 147 : dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
159 147 : if( !pbc[i] ) {
160 95 : max[i] +=dx[i];
161 95 : nbin[i]+=1;
162 : }
163 147 : stride[i]=npoints;
164 147 : npoints*=nbin[i];
165 : }
166 98 : resize(); // Always resize after setting new bounds as grid size may have have changed
167 98 : }
168 :
169 3 : void GridVessel::setupFibonacciGrid( const unsigned& np ) {
170 3 : bounds_set=true;
171 3 : root5 = std::sqrt(5);
172 3 : npoints = np;
173 3 : golden = ( 1 + std::sqrt(5) ) / 2.0;
174 3 : igolden = golden - 1;
175 3 : fib_increment = 2*pi*igolden;
176 3 : log_golden2 = std::log( golden*golden );
177 3 : fib_offset = 2 / static_cast<double>( npoints );
178 3 : fib_shift = fib_offset/2 - 1;
179 3 : resize();
180 :
181 3 : std::vector<double> icoord( dimension ), jcoord( dimension );
182 : // Find minimum distance between each pair of points
183 3 : std::vector<double> mindists( npoints );
184 347 : for(unsigned i=0; i<npoints; ++i) {
185 344 : getFibonacciCoordinates( i, icoord );
186 344 : mindists[i] = 0;
187 41080 : for(unsigned j=0; j<npoints; ++j) {
188 40736 : if( i==j ) {
189 344 : continue ; // Points are not neighbors to themselves
190 : }
191 40392 : getFibonacciCoordinates( j, jcoord );
192 : // Calculate the dot product
193 : double dot=0;
194 161568 : for(unsigned k=0; k<dimension; ++k) {
195 121176 : dot += icoord[k]*jcoord[k];
196 : }
197 40392 : if( dot>mindists[i] ) {
198 3055 : mindists[i]=dot;
199 : }
200 : }
201 : }
202 : // And now take minimum of dot products
203 3 : double min=mindists[0];
204 344 : for(unsigned i=1; i<npoints; ++i) {
205 341 : if( mindists[i]<min ) {
206 : min=mindists[i];
207 : }
208 : }
209 : double final_cutoff;
210 3 : if( getFibonacciCutoff()<-1 ) {
211 : final_cutoff=-1;
212 : } else {
213 2 : final_cutoff = std::cos( std::acos( getFibonacciCutoff() ) + std::acos( min ) );
214 : }
215 :
216 : // And now construct the neighbor list
217 3 : fib_nlist.resize( npoints );
218 347 : for(unsigned i=0; i<npoints; ++i) {
219 344 : getFibonacciCoordinates( i, icoord );
220 41080 : for(unsigned j=0; j<npoints; ++j) {
221 40736 : if( i==j ) {
222 344 : continue ; // Points are not neighbors to themselves
223 : }
224 40392 : getFibonacciCoordinates( j, jcoord );
225 : // Calculate the dot product
226 : double dot=0;
227 161568 : for(unsigned k=0; k<dimension; ++k) {
228 121176 : dot += icoord[k]*jcoord[k];
229 : }
230 40392 : if( dot>final_cutoff ) {
231 22886 : fib_nlist[i].push_back(j);
232 : }
233 : }
234 : }
235 3 : }
236 :
237 65 : std::string GridVessel::description() {
238 65 : if( !bounds_set ) {
239 14 : return "";
240 : }
241 :
242 : std::string des;
243 51 : if( gtype==flat ) {
244 : des="grid of ";
245 : std::string num;
246 68 : for(unsigned i=0; i<dimension-1; ++i) {
247 19 : Tools::convert( nbin[i], num );
248 38 : des += num + " X ";
249 : }
250 49 : Tools::convert( nbin[dimension-1], num );
251 49 : des += num + " equally spaced points between (";
252 68 : for(unsigned i=0; i<dimension-1; ++i) {
253 38 : des += str_min[i] + ",";
254 : }
255 49 : Tools::convert( nbin[dimension-1], num );
256 49 : des += str_min[dimension-1] + ") and (";
257 68 : for(unsigned i=0; i<dimension-1; ++i) {
258 38 : des += str_max[i] + ",";
259 : }
260 98 : des += str_max[dimension-1] + ")";
261 2 : } else if( gtype==fibonacci ) {
262 : std::string num;
263 2 : Tools::convert( npoints, num );
264 4 : des += "fibonacci grid of " + num + " points on spherical surface";
265 : }
266 : return des;
267 : }
268 :
269 213 : void GridVessel::resize() {
270 213 : plumed_massert( nper>0, "Number of datapoints at each grid point has not been set");
271 213 : if( getAction() ) {
272 210 : resizeBuffer( getNumberOfBufferPoints()*nper + 1 + 2*getAction()->getNumberOfDerivatives() );
273 : }
274 213 : setDataSize( npoints*nper );
275 213 : forces.resize( npoints );
276 213 : if( active.size()!=npoints) {
277 67 : active.resize( npoints, true );
278 : }
279 213 : }
280 :
281 24333366 : unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
282 : plumed_dbg_assert( gtype==flat && bounds_set && indices.size()==dimension );
283 : // indices are flattended using a column-major order
284 24333366 : unsigned index=indices[dimension-1];
285 69887384 : for(unsigned i=dimension-1; i>0; --i) {
286 45554018 : index=index*nbin[i-1]+indices[i-1];
287 : }
288 24333366 : return index;
289 : }
290 :
291 1060387 : void GridVessel::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
292 : plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension && indices.size()==dimension );
293 3196126 : for(unsigned i=0; i<dimension; ++i) {
294 2135739 : indices[i]=std::floor( (point[i] - min[i])/dx[i] );
295 2135739 : if( pbc[i] ) {
296 45516 : indices[i]=indices[i]%nbin[i];
297 2090223 : } else if( indices[i]>nbin[i] ) {
298 : std::string pp, num;
299 0 : Tools::convert( point[0], pp );
300 0 : for(unsigned j=1; j<point.size(); ++j) {
301 0 : Tools::convert( point[j], num );
302 0 : pp += ", " + num;
303 : }
304 0 : plumed_merror("point (" + pp + ") is outside grid range");
305 : }
306 : }
307 1060387 : }
308 :
309 530955 : unsigned GridVessel::getIndex( const std::vector<double>& point ) const {
310 : plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension );
311 530955 : if( gtype==flat ) {
312 530955 : std::vector<unsigned> indices(dimension);
313 530955 : getIndices( point, indices );
314 530955 : return getIndex( indices );
315 0 : } else if( gtype==fibonacci ) {
316 0 : return getFibonacciIndex( point );
317 : } else {
318 0 : plumed_error();
319 : }
320 : }
321 :
322 57 : unsigned GridVessel::getFibonacciIndex( const std::vector<double>& p ) const {
323 : plumed_dbg_assert( gtype==fibonacci );
324 : // Convert input point to coordinates on cylinder
325 : int k=2;
326 57 : double phi = std::atan2( p[2], p[0] ), sinthet2 = 1 - p[1]*p[1];
327 : // Calculate power to raise golden ratio
328 57 : if( sinthet2<epsilon ) {
329 : k = 2;
330 : } else {
331 57 : k = std::floor( std::log( npoints*pi*root5*sinthet2 ) / log_golden2 );
332 : if( k<2 ) {
333 : k = 2;
334 : }
335 : }
336 57 : double Fk = pow( golden, k ) / root5, F0 = std::round(Fk), F1 = std::round(Fk*golden);
337 : Matrix<double> B(2,2), invB(2,2);
338 57 : std::vector<double> thisp(3);
339 57 : B(0,0) = 2*pi*((F0+1)*igolden - std::floor((F0+1)*igolden)) - fib_increment;
340 57 : B(0,1) = 2*pi*((F1+1)*igolden - std::floor((F1+1)*igolden)) - fib_increment;
341 57 : B(1,0) = -2*F0/npoints;
342 57 : B(1,1) = -2*F1/npoints;
343 57 : Invert( B, invB );
344 57 : std::vector<double> vv(2), rc(2);
345 57 : vv[0]=-phi;
346 57 : vv[1] = p[1] - fib_shift;
347 57 : mult( invB, vv, rc );
348 57 : std::vector<int> c(2);
349 57 : c[0]=std::floor(rc[0]);
350 57 : c[1]=std::floor(rc[1]);
351 : unsigned outind=0;
352 : double mindist = 10000000.;
353 285 : for(int s=0; s<4; ++s) {
354 228 : double ttt, costheta = B(1,0)*( c[0] + s%2 ) + B(1,1)*( c[1] + s/2 ) + fib_shift;
355 228 : if( costheta>1 ) {
356 : ttt=1;
357 225 : } else if( costheta<-1 ) {
358 : ttt=-1;
359 : } else {
360 : ttt=costheta;
361 : }
362 228 : costheta = 2*ttt - costheta;
363 228 : unsigned i = std::floor( 0.5*npoints*(1+costheta) );
364 228 : getFibonacciCoordinates( i, thisp );
365 : double dist=0;
366 912 : for(unsigned j=0; j<3; ++j) {
367 684 : double tmp=thisp[j]-p[j];
368 684 : dist += tmp*tmp;
369 : }
370 228 : if( dist<mindist ) {
371 132 : outind = i;
372 : mindist = dist;
373 : }
374 : }
375 57 : return outind;
376 : }
377 :
378 112427807 : void GridVessel::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
379 : plumed_dbg_assert( gtype==flat );
380 112427807 : unsigned kk=index;
381 112427807 : indices[0]=index%nnbin[0];
382 220890741 : for(unsigned i=1; i<dimension-1; ++i) {
383 108462934 : kk=(kk-indices[i-1])/nnbin[i-1];
384 108462934 : indices[i]=kk%nnbin[i];
385 : }
386 112427807 : if(dimension>=2) { // I think this is wrong
387 112411300 : indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
388 : }
389 112427807 : }
390 :
391 15879009 : void GridVessel::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
392 : plumed_dbg_assert( gtype==flat );
393 15879009 : convertIndexToIndices( index, nbin, indices );
394 15879009 : }
395 :
396 665380 : void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
397 665380 : std::vector<unsigned> tindices( dimension );
398 665380 : getGridPointCoordinates( ipoint, tindices, x );
399 665380 : }
400 :
401 12514073 : void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
402 : plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
403 12514073 : if( gtype==flat ) {
404 12508956 : getFlatGridCoordinates( ipoint, tindices, x );
405 5117 : } else if( gtype==fibonacci ) {
406 5117 : getFibonacciCoordinates( ipoint, x );
407 : } else {
408 0 : plumed_error();
409 : }
410 12514073 : }
411 :
412 13037800 : void GridVessel::getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
413 : plumed_dbg_assert( gtype==flat );
414 13037800 : getIndices( ipoint, tindices );
415 50950441 : for(unsigned i=0; i<dimension; ++i) {
416 37912641 : x[i] = min[i] + dx[i]*tindices[i];
417 : }
418 13037800 : }
419 :
420 86817 : void GridVessel::getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
421 : plumed_dbg_assert( gtype==fibonacci );
422 86817 : x[1] = (ipoint*fib_offset) + fib_shift;
423 86817 : double r = std::sqrt( 1 - x[1]*x[1] );
424 86817 : double phi = ipoint*fib_increment;
425 86817 : x[0] = r*std::cos(phi);
426 86817 : x[2] = r*std::sin(phi);
427 : double norm=0;
428 347268 : for(unsigned j=0; j<3; ++j) {
429 260451 : norm+=x[j]*x[j];
430 : }
431 86817 : norm = std::sqrt(norm);
432 347268 : for(unsigned j=0; j<3; ++j) {
433 260451 : x[j] = x[j] / norm;
434 : }
435 86817 : }
436 :
437 528844 : void GridVessel::getSplineNeighbors( const unsigned& mybox, unsigned& nneighbors, std::vector<unsigned>& mysneigh ) const {
438 : plumed_dbg_assert( gtype==flat );
439 528844 : unsigned nneigh=unsigned(pow(2.0,int(dimension)));
440 528844 : if( mysneigh.size()!=nneigh ) {
441 528844 : mysneigh.resize(nneigh);
442 : }
443 :
444 : unsigned inind;
445 528844 : nneighbors = 0;
446 528844 : std::vector<unsigned> tmp_indices( dimension );
447 528844 : std::vector<unsigned> my_indices( dimension );
448 528844 : getIndices( mybox, my_indices );
449 2677596 : for(unsigned i=0; i<nneigh; ++i) {
450 : unsigned tmp=i;
451 : inind=0;
452 6513408 : for(unsigned j=0; j<dimension; ++j) {
453 4364656 : unsigned i0=tmp%2+my_indices[j];
454 4364656 : tmp/=2;
455 4364656 : if(!pbc[j] && i0==nbin[j]) {
456 40800 : continue;
457 : }
458 4323856 : if( pbc[j] && i0==nbin[j]) {
459 : i0=0;
460 : }
461 4323856 : tmp_indices[inind++]=i0;
462 : }
463 2148752 : if(inind==dimension ) {
464 2108152 : unsigned findex=getIndex( tmp_indices );
465 2108152 : mysneigh[nneighbors++]=findex;
466 2108152 : plumed_massert( active[findex], "inactive grid point required for splines");
467 : }
468 : }
469 528844 : }
470 :
471 6552681 : double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
472 13105362 : plumed_assert( bounds_set && ipoint<npoints && jelement<nper && active[ipoint] );
473 6552681 : return getDataElement( nper*ipoint + jelement );
474 : }
475 :
476 154208 : void GridVessel::setGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ) {
477 : plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
478 154208 : setDataElement( nper*ipoint + jelement, value );
479 154208 : }
480 :
481 24200 : void GridVessel::setValueAndDerivatives( const unsigned& ipoint, const unsigned& jelement, const double& value, const std::vector<double>& der ) {
482 : plumed_dbg_assert( !noderiv && jelement<getNumberOfComponents() && der.size()==nbin.size() );
483 24200 : setGridElement( ipoint, jelement, value );
484 72600 : for(unsigned i=0; i<der.size(); ++i) {
485 48400 : setGridElement( ipoint, jelement+1+i, der[i] );
486 : }
487 24200 : }
488 :
489 2605 : void GridVessel::addToGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ) {
490 : plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
491 2605 : addDataElement( nper*ipoint + jelement, value );
492 2605 : }
493 :
494 15087 : void GridVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
495 : plumed_dbg_assert( myvals.getNumberOfValues()==(nper+1) );
496 65138 : for(unsigned i=0; i<nper; ++i) {
497 50051 : buffer[bufstart + nper*current + i] += myvals.get(i+1);
498 : }
499 15087 : }
500 :
501 118 : void GridVessel::finish( const std::vector<double>& buffer ) {
502 118 : if( wasforced ) {
503 20 : getFinalForces( buffer, finalForces );
504 : } else {
505 98 : AveragingVessel::finish( buffer );
506 : }
507 118 : }
508 :
509 0 : double GridVessel::getGridElement( const std::vector<unsigned>& indices, const unsigned& jelement ) const {
510 0 : return getGridElement( getIndex( indices ), jelement );
511 : }
512 :
513 0 : void GridVessel::setGridElement( const std::vector<unsigned>& indices, const unsigned& jelement, const double& value ) {
514 0 : setGridElement( getIndex( indices ), jelement, value );
515 0 : }
516 :
517 202520 : std::vector<std::string> GridVessel::getMin() const {
518 : plumed_dbg_assert( gtype==flat );
519 202520 : return str_min;
520 : }
521 :
522 202553 : std::vector<std::string> GridVessel::getMax() const {
523 : plumed_dbg_assert( gtype==flat );
524 202553 : return str_max;
525 : }
526 :
527 203150 : std::vector<unsigned> GridVessel::getNbin() const {
528 : plumed_dbg_assert( gtype==flat && bounds_set );
529 203150 : std::vector<unsigned> ngrid( dimension );
530 605377 : for(unsigned i=0; i<dimension; ++i) {
531 402227 : if( !pbc[i] ) {
532 368027 : ngrid[i]=nbin[i] - 1;
533 : } else {
534 34200 : ngrid[i]=nbin[i];
535 : }
536 : }
537 203150 : return ngrid;
538 : }
539 :
540 21514 : void GridVessel::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
541 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
542 : plumed_dbg_assert( bounds_set );
543 :
544 21514 : if( gtype == flat ) {
545 : plumed_dbg_assert( nneigh.size()==dimension );
546 21457 : std::vector<unsigned> indices( dimension );
547 85340 : for(unsigned i=0; i<dimension; ++i) {
548 63883 : indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
549 : }
550 21457 : getNeighbors( indices, nneigh, num_neighbors, neighbors );
551 57 : } else if( gtype == fibonacci ) {
552 57 : unsigned find = getFibonacciIndex( pp );
553 57 : num_neighbors = 1 + fib_nlist[find].size();
554 57 : if( neighbors.size()<num_neighbors ) {
555 57 : neighbors.resize( num_neighbors );
556 : }
557 57 : neighbors[0]=find;
558 4773 : for(unsigned i=0; i<fib_nlist[find].size(); ++i) {
559 4716 : neighbors[1+i] = fib_nlist[find][i];
560 : }
561 : } else {
562 0 : plumed_error();
563 : }
564 21514 : }
565 :
566 40878 : void GridVessel::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
567 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
568 : plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
569 :
570 : unsigned num_neigh=1;
571 40878 : std::vector<unsigned> small_bin( dimension );
572 163024 : for(unsigned i=0; i<dimension; ++i) {
573 122146 : small_bin[i]=(2*nneigh[i]+1);
574 122146 : num_neigh *=small_bin[i];
575 : }
576 40878 : if( neighbors.size()!=num_neigh ) {
577 22050 : neighbors.resize( num_neigh );
578 : }
579 :
580 40878 : num_neighbors=0;
581 40878 : std::vector<unsigned> s_indices(dimension), t_indices(dimension);
582 96589676 : for(unsigned index=0; index<num_neigh; ++index) {
583 : bool found=true;
584 96548798 : convertIndexToIndices( index, small_bin, s_indices );
585 386166232 : for(unsigned i=0; i<dimension; ++i) {
586 289617434 : int i0=s_indices[i]-nneigh[i]+indices[i];
587 289617434 : if(!pbc[i] && i0<0) {
588 : found=false;
589 : }
590 289617434 : if(!pbc[i] && i0>=nbin[i]) {
591 : found=false;
592 : }
593 289617434 : if( pbc[i] && i0<0) {
594 15055490 : i0=nbin[i]-(-i0)%nbin[i];
595 : }
596 289617434 : if( pbc[i] && i0>=nbin[i]) {
597 15689040 : i0%=nbin[i];
598 : }
599 289617434 : t_indices[i]=static_cast<unsigned>(i0);
600 : }
601 96548798 : if( found ) {
602 21093803 : neighbors[num_neighbors]=getIndex( t_indices );
603 21093803 : num_neighbors++;
604 : }
605 : }
606 40878 : }
607 :
608 11 : void GridVessel::setCubeUnits( const double& units ) {
609 : plumed_dbg_assert( gtype==flat );
610 11 : cube_units=units;
611 11 : }
612 :
613 8 : double GridVessel::getCubeUnits() const {
614 : plumed_dbg_assert( gtype==flat );
615 8 : return cube_units;
616 : }
617 :
618 17 : std::string GridVessel::getInputString() const {
619 17 : std::string mstring="COORDINATES="+arg_names[0];
620 23 : for(unsigned i=1; i<dimension; ++i) {
621 12 : mstring+="," + arg_names[i];
622 : }
623 17 : if( gtype==flat ) {
624 : mstring += " TYPE=flat PBC=";
625 17 : if( pbc[0] ) {
626 : mstring +="T";
627 : } else {
628 : mstring +="F";
629 : }
630 23 : for(unsigned i=1; i<dimension; ++i) {
631 6 : if( pbc[i] ) {
632 : mstring +=",T";
633 : } else {
634 : mstring +=",F";
635 : }
636 : }
637 0 : } else if( gtype==fibonacci ) {
638 : mstring += " TYPE=fibonacci";
639 : }
640 17 : return mstring;
641 : }
642 :
643 528844 : double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const {
644 : plumed_dbg_assert( gtype==flat && der.size()==dimension && !noderiv && ind<getNumberOfComponents() );
645 :
646 : double X,X2,X3,value=0;
647 528844 : der.assign(der.size(),0.0);
648 528844 : std::vector<double> fd(dimension);
649 528844 : std::vector<double> C(dimension);
650 528844 : std::vector<double> D(dimension);
651 528844 : std::vector<double> dder(dimension);
652 :
653 528844 : std::vector<unsigned> nindices(dimension);
654 : unsigned n_neigh;
655 528844 : std::vector<unsigned> indices(dimension);
656 528844 : getIndices( x, indices );
657 : std::vector<unsigned> neigh;
658 528844 : getSplineNeighbors( getIndex(indices), n_neigh, neigh );
659 528844 : std::vector<double> xfloor(dimension);
660 528844 : getFlatGridCoordinates( getIndex(x), nindices, xfloor );
661 :
662 : // loop over neighbors
663 2636996 : for(unsigned int ipoint=0; ipoint<n_neigh; ++ipoint) {
664 2108152 : double grid=getGridElement(neigh[ipoint], ind*(1+dimension) );
665 6391608 : for(unsigned j=0; j<dimension; ++j) {
666 4283456 : dder[j] = getGridElement( neigh[ipoint], ind*(1+dimension) + 1 + j );
667 : }
668 :
669 2108152 : getIndices( neigh[ipoint], nindices );
670 : double ff=1.0;
671 :
672 6391608 : for(unsigned j=0; j<dimension; ++j) {
673 : int x0=1;
674 4283456 : if(nindices[j]==indices[j]) {
675 : x0=0;
676 : }
677 4283456 : double ddx=dx[j];
678 4283456 : X=std::fabs((x[j]-xfloor[j])/ddx-(double)x0);
679 4283456 : X2=X*X;
680 4283456 : X3=X2*X;
681 : double yy;
682 4283456 : if(std::fabs(grid)<0.0000001) {
683 : yy=0.0;
684 : } else {
685 4269609 : yy=-dder[j]/grid;
686 : }
687 6445384 : C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*ddx;
688 4283456 : D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*ddx;
689 4283456 : D[j]*=(x0?-1.0:1.0)/ddx;
690 4283456 : ff*=C[j];
691 : }
692 6391608 : for(unsigned j=0; j<dimension; ++j) {
693 4283456 : fd[j]=D[j];
694 13052624 : for(unsigned i=0; i<dimension; ++i)
695 8769168 : if(i!=j) {
696 4485712 : fd[j]*=C[i];
697 : }
698 : }
699 2108152 : value+=grid*ff;
700 6391608 : for(unsigned j=0; j<dimension; ++j) {
701 4283456 : der[j]+=grid*fd[j];
702 : }
703 : }
704 528844 : return value;
705 : }
706 :
707 3 : void GridVessel::activateThesePoints( const std::vector<bool>& to_activate ) {
708 : plumed_dbg_assert( to_activate.size()==npoints );
709 29991 : for(unsigned i=0; i<npoints; ++i) {
710 29988 : active[i]=to_activate[i];
711 : }
712 3 : }
713 :
714 20 : void GridVessel::setForce( const std::vector<double>& inforces ) {
715 : plumed_dbg_assert( inforces.size()==npoints );
716 20 : wasforced=true;
717 5725 : for(unsigned i=0; i<npoints; ++i) {
718 5705 : forces[i]=inforces[i];
719 : }
720 20 : }
721 :
722 11870273 : bool GridVessel::wasForced() const {
723 11870273 : return wasforced;
724 : }
725 :
726 20 : bool GridVessel::applyForce( std::vector<double>& fforces ) {
727 : plumed_dbg_assert( fforces.size()==finalForces.size() );
728 20 : if( !wasforced ) {
729 : return false;
730 : }
731 3815 : for(unsigned i=0; i<finalForces.size(); ++i) {
732 3795 : fforces[i]=finalForces[i];
733 : }
734 20 : wasforced=false;
735 20 : return true;
736 : }
737 :
738 : }
739 : }
740 :
|
