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 "GridCoordinatesObject.h"
23 : #include "tools/Random.h"
24 : #include "tools/Matrix.h"
25 :
26 : namespace PLMD {
27 : namespace gridtools {
28 :
29 335 : void GridCoordinatesObject::setup( const std::string& geom, const std::vector<bool>& ipbc,
30 : const unsigned& np, const double& fib_cutoff ) {
31 335 : if( geom=="flat" ) {
32 318 : gtype=flat;
33 318 : dimension = ipbc.size();
34 17 : } else if( geom=="fibonacci" ) {
35 17 : gtype=fibonacci;
36 17 : dimension = 3;
37 : } else {
38 0 : plumed_merror( geom + " is invalid geometry type");
39 : }
40 :
41 335 : if( gtype==flat ) {
42 318 : bounds_set=false;
43 318 : npoints=0;
44 318 : pbc.resize( ipbc.size() );
45 678 : for(unsigned i=0; i<ipbc.size(); ++i) {
46 360 : pbc[i]=ipbc[i];
47 : }
48 17 : } else if( gtype==fibonacci ) {
49 17 : bounds_set=true;
50 17 : root5 = sqrt(5);
51 17 : npoints = np;
52 17 : golden = ( 1 + sqrt(5) ) / 2.0;
53 17 : igolden = golden - 1;
54 17 : fib_increment = 2*pi*igolden;
55 17 : log_golden2 = std::log( golden*golden );
56 17 : fib_offset = 2 / static_cast<double>( npoints );
57 17 : fib_shift = fib_offset/2 - 1;
58 :
59 17 : std::vector<double> icoord( dimension ), jcoord( dimension );
60 : // Find minimum distance between each pair of points
61 17 : std::vector<unsigned> tindices( dimension );
62 17 : std::vector<double> mindists( npoints );
63 4849 : for(unsigned i=0; i<npoints; ++i) {
64 4832 : getFibonacciCoordinates( i, icoord );
65 4832 : mindists[i] = 0;
66 1707040 : for(unsigned j=0; j<npoints; ++j) {
67 1702208 : if( i==j ) {
68 4832 : continue ; // Points are not neighbors to themselves
69 : }
70 1697376 : getFibonacciCoordinates( j, jcoord );
71 : // Calculate the dot product
72 : double dot=0;
73 6789504 : for(unsigned k=0; k<dimension; ++k) {
74 5092128 : dot += icoord[k]*jcoord[k];
75 : }
76 1697376 : if( dot>mindists[i] ) {
77 85497 : mindists[i]=dot;
78 : }
79 : }
80 : }
81 : // And now take minimum of dot products
82 17 : double min=mindists[0];
83 4832 : for(unsigned i=1; i<npoints; ++i) {
84 4815 : if( mindists[i]<min ) {
85 : min=mindists[i];
86 : }
87 : }
88 : double final_cutoff;
89 17 : if( fib_cutoff<-1 ) {
90 : final_cutoff=-1;
91 : } else {
92 15 : final_cutoff = cos( acos( fib_cutoff ) + acos( min ) );
93 : }
94 :
95 : // And now construct the neighbor list
96 17 : fib_nlist.resize( npoints );
97 4849 : for(unsigned i=0; i<npoints; ++i) {
98 4832 : getFibonacciCoordinates( i, icoord );
99 1707040 : for(unsigned j=0; j<npoints; ++j) {
100 1702208 : if( i==j ) {
101 4832 : continue ; // Points are not neighbors to themselves
102 : }
103 1697376 : getFibonacciCoordinates( j, jcoord );
104 : // Calculate the dot product
105 : double dot=0;
106 6789504 : for(unsigned k=0; k<dimension; ++k) {
107 5092128 : dot += icoord[k]*jcoord[k];
108 : }
109 1697376 : if( dot>final_cutoff ) {
110 794204 : fib_nlist[i].push_back(j);
111 : }
112 : }
113 : }
114 : } else {
115 0 : plumed_error();
116 : }
117 335 : }
118 :
119 618 : void GridCoordinatesObject::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
120 : const std::vector<unsigned>& binsin, std::vector<double>& spacing ) {
121 : plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
122 618 : plumed_assert( gtype==flat && (spacing.size()==dimension || binsin.size()==dimension) );
123 618 : str_min.resize( dimension );
124 618 : str_max.resize( dimension );
125 618 : nbin.resize( dimension );
126 618 : min.resize( dimension );
127 618 : max.resize( dimension );
128 618 : dx.resize( dimension );
129 618 : stride.resize( dimension );
130 :
131 618 : npoints=1;
132 1164 : bounds_set=(smin[0]!="auto" && smax[0]!="auto");
133 618 : if( bounds_set ) {
134 596 : for(unsigned i=1; i<dimension; ++i) {
135 100 : if( smin[i]=="auto" || smax[i]=="auto" ) {
136 0 : bounds_set=false;
137 0 : break;
138 : }
139 : }
140 : }
141 1302 : for(unsigned i=0; i<dimension; ++i) {
142 684 : str_min[i]=smin[i];
143 : str_max[i]=smax[i];
144 684 : if( bounds_set ) {
145 596 : Tools::convert( str_min[i], min[i] );
146 596 : Tools::convert( str_max[i], max[i] );
147 : }
148 684 : if( spacing.size()==dimension && binsin.size()==dimension ) {
149 290 : if( spacing[i]==0 ) {
150 17 : nbin[i] = binsin[i];
151 273 : } else if( bounds_set ) {
152 273 : double range = max[i] - min[i];
153 273 : nbin[i] = std::round( range / spacing[i]);
154 273 : dx[i]=spacing[i];
155 : // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
156 273 : if( nbin[i]!=binsin[i] ) {
157 0 : plumed_merror("mismatch between input spacing and input number of bins");
158 : }
159 : }
160 394 : } else if( binsin.size()==dimension ) {
161 394 : nbin[i]=binsin[i];
162 394 : dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
163 0 : } else if( spacing.size()==dimension && bounds_set ) {
164 0 : nbin[i] = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
165 0 : dx[i]=spacing[i];
166 0 : } else if( bounds_set ) {
167 0 : plumed_error();
168 : }
169 684 : if( !pbc[i] ) {
170 481 : max[i] +=dx[i];
171 481 : nbin[i]+=1;
172 : }
173 684 : stride[i]=npoints;
174 684 : npoints*=nbin[i];
175 : }
176 618 : if( spacing.size()!=dimension && bounds_set ) {
177 293 : spacing.resize(dimension);
178 616 : for(unsigned i=0; i<dimension; ++i) {
179 323 : spacing[i]=dx[i];
180 : }
181 : }
182 618 : }
183 :
184 42783311 : unsigned GridCoordinatesObject::getIndex( const std::vector<unsigned>& indices ) const {
185 : plumed_dbg_assert( gtype==flat && indices.size()==dimension );
186 : // indices are flattended using a column-major order
187 42783311 : unsigned index=indices[dimension-1];
188 122880687 : for(unsigned i=dimension-1; i>0; --i) {
189 80097376 : index=index*nbin[i-1]+indices[i-1];
190 : }
191 42783311 : return index;
192 : }
193 :
194 197426 : bool GridCoordinatesObject::inbounds( const std::vector<double>& point ) const {
195 197426 : if( gtype==fibonacci ) {
196 : return true;
197 : }
198 : plumed_dbg_assert( bounds_set && point.size()==dimension );
199 461467 : for(unsigned i=0; i<dimension; ++i) {
200 296963 : if( pbc[i] ) {
201 236124 : continue;
202 : }
203 60839 : if( point[i]<min[i] || point[i]>(max[i]-dx[i]) ) {
204 : return false;
205 : }
206 : }
207 : return true;
208 : }
209 :
210 1175406 : void GridCoordinatesObject::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
211 : plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension && indices.size()==dimension );
212 3445432 : for(unsigned i=0; i<dimension; ++i) {
213 2270026 : indices[i]=std::floor( (point[i] - min[i])/dx[i] );
214 2270026 : if( pbc[i] ) {
215 90308 : indices[i]=indices[i]%nbin[i];
216 2179718 : } else if( indices[i]>nbin[i] ) {
217 0 : plumed_merror("point is outside grid range");
218 : }
219 : }
220 1175406 : }
221 :
222 615967 : unsigned GridCoordinatesObject::getIndex( const std::vector<double>& point ) const {
223 : plumed_dbg_assert( bounds_set && point.size()==dimension );
224 615967 : if( gtype==flat ) {
225 613867 : std::vector<unsigned> indices(dimension);
226 613867 : getIndices( point, indices );
227 613867 : return getIndex( indices );
228 2100 : } else if( gtype==fibonacci ) {
229 2100 : return getFibonacciIndex( point );
230 : } else {
231 0 : plumed_error();
232 : }
233 : }
234 :
235 16510 : unsigned GridCoordinatesObject::getFibonacciIndex( const std::vector<double>& p ) const {
236 : plumed_dbg_assert( gtype==fibonacci );
237 : // Convert input point to coordinates on cylinder
238 : int k=2;
239 16510 : double phi = atan2( p[2], p[0] ), sinthet2 = 1 - p[1]*p[1];
240 : // Calculate power to raise golden ratio
241 16510 : if( sinthet2<epsilon ) {
242 : k = 2;
243 : } else {
244 16510 : k = std::floor( std::log( npoints*pi*root5*sinthet2 ) / log_golden2 );
245 : if( k<2 ) {
246 : k = 2;
247 : }
248 : }
249 16510 : double Fk = pow( golden, k ) / root5, F0 = std::round(Fk), F1 = std::round(Fk*golden);
250 : Matrix<double> B(2,2), invB(2,2);
251 16510 : std::vector<double> thisp(3);
252 16510 : B(0,0) = 2*pi*((F0+1)*igolden - std::floor((F0+1)*igolden)) - fib_increment;
253 16510 : B(0,1) = 2*pi*((F1+1)*igolden - std::floor((F1+1)*igolden)) - fib_increment;
254 16510 : B(1,0) = -2*F0/npoints;
255 16510 : B(1,1) = -2*F1/npoints;
256 16510 : Invert( B, invB );
257 16510 : std::vector<double> vv(2), rc(2);
258 16510 : vv[0]=-phi;
259 16510 : vv[1] = p[1] - fib_shift;
260 16510 : mult( invB, vv, rc );
261 16510 : std::vector<int> c(2);
262 16510 : c[0]=std::floor(rc[0]);
263 16510 : c[1]=std::floor(rc[1]);
264 : unsigned outind=0;
265 : double mindist = 10000000.;
266 82550 : for(int s=0; s<4; ++s) {
267 66040 : double ttt, costheta = B(1,0)*( c[0] + s%2 ) + B(1,1)*( c[1] + s/2 ) + fib_shift;
268 66040 : if( costheta>1 ) {
269 : ttt=1;
270 64959 : } else if( costheta<-1 ) {
271 : ttt=-1;
272 : } else {
273 : ttt=costheta;
274 : }
275 66040 : costheta = 2*ttt - costheta;
276 66040 : unsigned i = std::floor( 0.5*npoints*(1+costheta) );
277 66040 : getFibonacciCoordinates( i, thisp );
278 : double dist=0;
279 264160 : for(unsigned j=0; j<3; ++j) {
280 198120 : double tmp=thisp[j]-p[j];
281 198120 : dist += tmp*tmp;
282 : }
283 66040 : if( dist<mindist ) {
284 33971 : outind = i;
285 : mindist = dist;
286 : }
287 : }
288 16510 : return outind;
289 : }
290 :
291 94917677 : void GridCoordinatesObject::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
292 : plumed_dbg_assert( gtype==flat );
293 94917677 : unsigned kk=index;
294 94917677 : indices[0]=index%nnbin[0];
295 183492225 : for(unsigned i=1; i<dimension-1; ++i) {
296 88574548 : kk=(kk-indices[i-1])/nnbin[i-1];
297 88574548 : indices[i]=kk%nnbin[i];
298 : }
299 94917677 : if(dimension>=2) { // I think this is wrong
300 92728227 : indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
301 : }
302 94917677 : }
303 :
304 42760057 : void GridCoordinatesObject::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
305 : plumed_dbg_assert( gtype==flat );
306 42760057 : convertIndexToIndices( index, nbin, indices );
307 42760057 : }
308 :
309 40720264 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
310 40720264 : std::vector<unsigned> tindices( dimension );
311 40720264 : getGridPointCoordinates( ipoint, tindices, x );
312 40720264 : }
313 :
314 41484313 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
315 : plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
316 41484313 : if( gtype==flat ) {
317 39980732 : getFlatGridCoordinates( ipoint, tindices, x );
318 1503581 : } else if( gtype==fibonacci ) {
319 1503581 : getFibonacciCoordinates( ipoint, x );
320 : } else {
321 0 : plumed_error();
322 : }
323 41484313 : }
324 :
325 0 : void GridCoordinatesObject::putCoordinateAtValue( const unsigned& ind, const double& val, std::vector<double>& coords ) const {
326 0 : std::vector<double> point( dimension );
327 0 : getGridPointCoordinates( ind, point );
328 0 : if( gtype==flat ) {
329 0 : if( coords.size()!=(dimension+1) ) {
330 0 : coords.resize( (dimension+1) );
331 : }
332 0 : for(unsigned i=0; i<dimension; ++i) {
333 0 : coords[i]=point[i];
334 : }
335 0 : coords[point.size()]=val;
336 0 : } else if( gtype==fibonacci ) {
337 0 : if( coords.size()!=3 ) {
338 0 : coords.resize(3);
339 : }
340 0 : for(unsigned i=0; i<3; ++i) {
341 0 : coords[i] = val*point[i];
342 : }
343 : } else {
344 0 : plumed_error();
345 : }
346 0 : }
347 :
348 39980732 : void GridCoordinatesObject::getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
349 : plumed_dbg_assert( gtype==flat );
350 39980732 : getIndices( ipoint, tindices );
351 156503342 : for(unsigned i=0; i<dimension; ++i) {
352 116522610 : x[i] = min[i] + dx[i]*tindices[i];
353 : }
354 39980732 : }
355 :
356 4974037 : void GridCoordinatesObject::getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
357 : plumed_dbg_assert( gtype==fibonacci );
358 4974037 : x[1] = (ipoint*fib_offset) + fib_shift;
359 4974037 : double r = sqrt( 1 - x[1]*x[1] );
360 4974037 : double phi = ipoint*fib_increment;
361 4974037 : x[0] = r*cos(phi);
362 4974037 : x[2] = r*sin(phi);
363 : double norm=0;
364 19896148 : for(unsigned j=0; j<3; ++j) {
365 14922111 : norm+=x[j]*x[j];
366 : }
367 4974037 : norm = sqrt(norm);
368 19896148 : for(unsigned j=0; j<3; ++j) {
369 14922111 : x[j] = x[j] / norm;
370 : }
371 4974037 : }
372 :
373 544167 : void GridCoordinatesObject::getSplineNeighbors( const unsigned& mybox, unsigned& nneighbors, std::vector<unsigned>& mysneigh ) const {
374 : plumed_dbg_assert( gtype==flat );
375 544167 : mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
376 :
377 : unsigned inind;
378 544167 : nneighbors = 0;
379 544167 : std::vector<unsigned> tmp_indices( dimension );
380 544167 : std::vector<unsigned> my_indices( dimension );
381 544167 : getIndices( mybox, my_indices );
382 2754079 : for(unsigned i=0; i<mysneigh.size(); ++i) {
383 : unsigned tmp=i;
384 : inind=0;
385 6716896 : for(unsigned j=0; j<dimension; ++j) {
386 4506984 : unsigned i0=tmp%2+my_indices[j];
387 4506984 : tmp/=2;
388 4506984 : if(!pbc[j] && i0==nbin[j]) {
389 40800 : continue;
390 : }
391 4466184 : if( pbc[j] && i0==nbin[j]) {
392 : i0=0;
393 : }
394 4466184 : tmp_indices[inind++]=i0;
395 : }
396 2209912 : if( inind==dimension ) {
397 2169312 : mysneigh[nneighbors++]=getIndex( tmp_indices );
398 : }
399 : }
400 544167 : }
401 :
402 470843 : std::vector<std::string> GridCoordinatesObject::getMin() const {
403 : plumed_dbg_assert( gtype==flat );
404 470843 : return str_min;
405 : }
406 :
407 470843 : std::vector<std::string> GridCoordinatesObject::getMax() const {
408 : plumed_dbg_assert( gtype==flat );
409 470843 : return str_max;
410 : }
411 :
412 483329 : std::vector<unsigned> GridCoordinatesObject::getNbin( const bool& shape ) const {
413 : plumed_dbg_assert( gtype==flat && nbin.size()==dimension );
414 483329 : std::vector<unsigned> ngrid( dimension );
415 1591468 : for(unsigned i=0; i<dimension; ++i) {
416 1108139 : if( !pbc[i] && !shape ) {
417 574082 : ngrid[i]=nbin[i] - 1;
418 : } else {
419 534057 : ngrid[i]=nbin[i];
420 : }
421 : }
422 483329 : return ngrid;
423 : }
424 :
425 231602 : void GridCoordinatesObject::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
426 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
427 : plumed_dbg_assert( bounds_set );
428 :
429 231602 : if( gtype == flat ) {
430 : plumed_dbg_assert( nneigh.size()==dimension );
431 217192 : std::vector<unsigned> indices( dimension );
432 490941 : for(unsigned i=0; i<dimension; ++i) {
433 273749 : indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
434 : }
435 217192 : getNeighbors( indices, nneigh, num_neighbors, neighbors );
436 14410 : } else if( gtype == fibonacci ) {
437 14410 : unsigned find = getFibonacciIndex( pp );
438 14410 : num_neighbors = 1 + fib_nlist[find].size();
439 14410 : if( neighbors.size()<num_neighbors ) {
440 14410 : neighbors.resize( num_neighbors );
441 : }
442 14410 : neighbors[0]=find;
443 1500997 : for(unsigned i=0; i<fib_nlist[find].size(); ++i) {
444 1486587 : neighbors[1+i] = fib_nlist[find][i];
445 : }
446 : } else {
447 0 : plumed_error();
448 : }
449 231602 : }
450 :
451 245010 : void GridCoordinatesObject::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
452 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
453 : plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
454 :
455 : unsigned num_neigh=1;
456 245010 : std::vector<unsigned> small_bin( dimension );
457 602213 : for(unsigned i=0; i<dimension; ++i) {
458 357203 : small_bin[i]=(2*nneigh[i]+1);
459 357203 : if( pbc[i] && small_bin[i]>nbin[i] ) {
460 133400 : small_bin[i]=nbin[i];
461 : }
462 357203 : num_neigh *=small_bin[i];
463 : }
464 245010 : if( neighbors.size()!=num_neigh ) {
465 217689 : neighbors.resize( num_neigh );
466 : }
467 :
468 245010 : num_neighbors=0;
469 245010 : std::vector<unsigned> s_indices(dimension), t_indices(dimension);
470 52402630 : for(unsigned index=0; index<num_neigh; ++index) {
471 : bool found=true;
472 52157620 : convertIndexToIndices( index, small_bin, s_indices );
473 206161849 : for(unsigned i=0; i<dimension; ++i) {
474 154004229 : int i0=s_indices[i]-nneigh[i]+indices[i];
475 154004229 : if(!pbc[i] && i0<0) {
476 : found=false;
477 : }
478 154004229 : if(!pbc[i] && i0>=nbin[i]) {
479 : found=false;
480 : }
481 154004229 : if( pbc[i] && i0<0) {
482 8443398 : i0=nbin[i]-(-i0)%nbin[i];
483 : }
484 154004229 : if( pbc[i] && i0>=nbin[i]) {
485 8670218 : i0%=nbin[i];
486 : }
487 154004229 : t_indices[i]=static_cast<unsigned>(i0);
488 : }
489 52157620 : if( found ) {
490 39399065 : neighbors[num_neighbors]=getIndex( t_indices );
491 39399065 : num_neighbors++;
492 : }
493 : }
494 245010 : }
495 :
496 : }
497 : }
498 :
|
