From fd27f122df2fa5375680d8e13e375b544fea3bbe Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Fri, 5 May 2017 15:47:17 +0100
Subject: [PATCH] Added routines for truncating kernel evaluation on spherical
fibonacci grids
---
.../rt-dfg-wcsurf/mysurface.xyz.reference | 200 +++++++++---------
.../rt-spherical-integral/colvar.reference | 10 +-
.../rt-spherical-integral/forces.reference | 58 ++---
src/gridtools/GridVessel.cpp | 128 +++++++++--
src/gridtools/GridVessel.h | 19 +-
src/gridtools/HistogramOnGrid.cpp | 11 +-
src/gridtools/HistogramOnGrid.h | 1 +
7 files changed, 266 insertions(+), 161 deletions(-)
diff --git a/regtest/adjmat/rt-dfg-wcsurf/mysurface.xyz.reference b/regtest/adjmat/rt-dfg-wcsurf/mysurface.xyz.reference
index d68d0b82f..44746f761 100644
--- a/regtest/adjmat/rt-dfg-wcsurf/mysurface.xyz.reference
+++ b/regtest/adjmat/rt-dfg-wcsurf/mysurface.xyz.reference
@@ -1,102 +1,102 @@
100
Grid converted to xyz file
-X -20.8562 -274.8130 -33.1425
-X 66.4253 -274.8131 18.2044
-X -80.0064 -273.5955 41.0579
-X 37.7558 -274.8134 -101.8397
-X 47.2090 -274.8134 115.9678
-X -126.4082 -273.0455 -59.9072
-X 142.9675 -263.9868 -44.0784
-X -82.9297 -264.6495 141.5048
-X -34.1637 -242.1433 -159.0945
-X 142.4703 -241.3816 101.2043
-X -206.6447 -268.3255 25.7537
-X 141.5894 -263.6224 -166.3450
-X 6.7746 -210.7854 185.7717
-X -136.7488 -208.0931 -138.7649
-X 213.1740 -215.2100 10.8997
-X -171.2314 -211.7436 141.4800
-X 33.6909 -220.0528 -241.4799
-X 118.0120 -177.9869 171.3899
-X -239.7151 -199.5465 -55.1579
-X 212.5498 -188.0594 -120.4222
-X -80.8322 -191.4166 249.1663
-X -104.0823 -173.9927 -228.1912
-X 239.6331 -171.4259 101.6646
-X -241.6610 -160.2321 85.5914
-X 128.6707 -162.4754 -241.9468
-X 64.6759 -145.3199 250.3071
-X -224.2947 -142.4388 -145.7711
-X 241.4801 -123.3645 -40.2888
-X -175.4468 -135.6189 224.2733
-X -14.5224 -84.0142 -186.3329
-X 185.8863 -107.7394 173.6529
-X -274.8126 -109.4532 -2.6813
-X 209.6625 -105.3020 -188.3390
-X -26.8224 -96.5266 274.8133
-X -181.4369 -98.4902 -241.4959
-X 241.4801 -74.4417 45.1484
-X -243.1214 -80.3046 151.3425
-X 44.7517 -42.9526 -160.2221
-X 110.9327 -57.9740 218.7866
-X -274.8101 -63.0670 -103.4100
-X 241.4801 -50.3570 -96.9330
-X -131.8832 -52.5844 274.8101
-X -82.5840 -43.2106 -272.5754
-X 241.4801 -36.8033 143.1079
-X -274.8134 -31.0752 57.6120
-X 173.3147 -26.8605 -241.4801
-X 32.8887 -19.4218 274.8126
-X -241.4801 -15.9419 -207.5827
-X 241.4801 -7.2513 -7.6268
-X -241.0937 -3.3689 235.2836
-X 10.7868 1.9285 -192.5430
-X 170.2794 8.0714 208.1467
-X -274.8110 13.9006 -39.7089
-X 241.4509 20.5025 -164.5226
-X -64.6454 25.5112 274.8061
-X -153.9322 34.8602 -274.8134
-X 241.4800 33.3386 79.6427
-X -274.7596 45.7159 123.7074
-X 118.1954 51.6147 -274.8606
-X 241.4767 46.7318 0.0000
-X -242.7131 70.7003 222.3451
-X 24.0852 65.1089 -274.4382
-X 159.5822 67.7208 208.1467
-X -274.8134 78.2580 -48.6106
-X 241.4197 86.7022 -153.5714
-X -73.8736 92.7852 274.8060
-X -139.0544 105.4684 -267.7409
-X 241.4776 96.0521 88.1872
-X -269.7450 116.2225 111.3468
-X 114.1679 114.0849 -243.9701
-X 75.2354 113.0023 239.8624
-X -241.4818 132.9307 -139.9437
-X 241.4197 124.5574 -53.0754
-X -153.5999 142.2089 218.4801
-X -32.2901 141.2368 -249.1804
-X 205.1457 159.0684 172.8970
-X -244.0258 152.6469 10.0912
-X 188.5883 175.2119 -187.6705
-X -9.8697 148.2284 213.4405
-X -149.1039 170.0551 -178.6764
-X 225.8814 175.4530 30.3920
-X -182.3928 180.2542 126.9042
-X 47.4924 185.0818 -211.1084
-X 114.6629 208.1460 200.1021
-X -205.6525 205.7813 -65.6088
-X 179.1695 198.9938 -82.7809
-X -75.6512 209.3042 180.7648
-X -64.4692 215.9868 -179.2409
-X 173.8660 237.0384 91.3791
-X -178.4690 237.8158 47.0439
-X 83.1987 212.4801 -129.3930
-X 27.4856 241.4801 159.9310
-X -118.3214 241.4801 -91.6348
-X 136.3858 241.4801 -11.2993
-X -84.1522 241.8844 90.9661
-X 0.5523 248.4932 -113.2155
-X 72.9756 274.8134 80.4450
-X -80.9491 247.3379 -7.5023
-X 54.8628 274.8134 -41.6389
-X -7.0088 274.8134 38.5265
+X 39.1587 -274.8129 0.0000
+X -50.7787 -274.7733 -46.5174
+X 7.8969 -274.8134 89.9809
+X 66.0845 -274.8134 -86.1955
+X -122.5208 -273.0884 21.6722
+X 109.9911 -254.4508 69.9673
+X -37.9272 -257.7900 -141.0874
+X -73.8311 -258.4720 142.1573
+X 171.4871 -271.6723 -62.6268
+X -172.1950 -257.3091 -71.0796
+X 79.2668 -240.9763 169.3886
+X 65.3260 -263.4165 -208.2696
+X -192.6446 -252.4683 111.6415
+X 200.4365 -219.2019 44.0654
+X -118.6361 -207.9757 -168.7477
+X -28.0598 -208.1467 216.5352
+X 177.6362 -209.6665 -149.7120
+X -239.9099 -205.3792 -9.9210
+X 147.6486 -168.9793 146.9301
+X -11.1467 -185.7668 -241.0569
+X -154.7927 -176.5432 185.4933
+X 241.4801 -169.0316 -32.4908
+X -208.1398 -166.9850 -144.8182
+X 55.0272 -156.6972 244.6015
+X 138.3903 -165.0346 -241.5095
+X -245.6773 -144.9540 78.3778
+X 234.0044 -137.2586 108.1160
+X -107.0039 -139.6660 -255.6804
+X -98.8447 -139.0974 274.8134
+X 241.4801 -122.6291 -126.9152
+X -274.8101 -120.3683 -72.4391
+X 133.8367 -98.5551 208.1467
+X 27.7848 -61.2915 -161.6720
+X -226.4329 -100.1200 175.3624
+X 241.4801 -79.0075 20.0061
+X -208.8177 -93.1796 -225.7258
+X 1.3406 -77.0624 274.8127
+X 200.0394 -76.8732 -220.5143
+X -274.8133 -65.2266 25.4696
+X 230.3384 -62.1100 174.8186
+X -36.7649 -39.7519 -202.0922
+X -155.7205 -50.4379 247.4554
+X 241.4801 -37.9875 -66.1797
+X -274.8105 -40.4987 -141.0282
+X 89.9687 -28.6436 242.6752
+X 71.6911 -17.1824 -176.1075
+X -274.8016 -21.3393 130.2336
+X 241.4801 -12.6507 74.4509
+X -161.0559 -9.5602 -274.8134
+X -59.0121 -2.8109 274.8089
+X 241.4791 2.9622 -171.5354
+X -274.8108 8.3118 -34.2491
+X 177.1701 13.6841 208.1467
+X 7.4913 14.4246 -205.4239
+X -236.9471 30.5054 240.4404
+X 241.4810 26.7600 -12.3471
+X -241.4828 41.0707 -199.5251
+X 38.3043 42.0567 274.5465
+X 166.2744 50.5786 -241.4819
+X -274.8134 54.5731 63.2339
+X 241.4801 59.6134 136.8130
+X -89.0838 68.2280 -274.6021
+X -125.3442 77.9870 274.8060
+X 241.4197 73.5381 -102.4226
+X -274.8134 88.3433 -97.3333
+X 113.9749 79.1471 214.3135
+X 70.9584 99.1560 -274.6213
+X -250.5401 111.6426 162.8282
+X 241.4775 97.5011 40.2884
+X -167.0284 114.8139 -213.5121
+X -19.4419 112.4741 249.4533
+X 218.6379 142.5028 -204.2490
+X -273.4190 137.7833 2.6677
+X 194.9349 139.5283 175.1093
+X -24.4079 141.2368 -250.0756
+X -148.5398 146.6198 197.7093
+X 241.4197 153.5023 -45.1372
+X -223.1565 173.1080 -138.9144
+X 58.9344 151.9795 210.9997
+X 109.8473 177.4985 -216.6460
+X -217.0556 178.5305 81.6772
+X 204.8306 179.0529 82.2214
+X -90.1106 178.1411 -187.7671
+X -63.0636 196.2906 208.1466
+X 175.2033 194.1460 -103.8304
+X -206.1657 212.3819 -43.2207
+X 125.5594 230.0046 174.9424
+X 21.8886 208.8659 -182.8973
+X -133.9740 213.2089 115.1677
+X 175.4127 226.1354 5.5401
+X -116.4589 224.7626 -113.6524
+X 9.0769 241.4801 162.0215
+X 89.7391 228.6845 -109.6956
+X -135.4464 241.4801 19.5713
+X 109.0941 257.6786 74.3358
+X -25.2617 242.1312 -107.3869
+X -48.4661 250.9330 86.5259
+X 85.7817 274.8134 -28.2917
+X -62.8027 274.8134 -28.2762
+X 15.4694 274.8134 35.9737
diff --git a/regtest/crystallization/rt-spherical-integral/colvar.reference b/regtest/crystallization/rt-spherical-integral/colvar.reference
index f8684e6a3..176f19c6f 100644
--- a/regtest/crystallization/rt-spherical-integral/colvar.reference
+++ b/regtest/crystallization/rt-spherical-integral/colvar.reference
@@ -1,6 +1,6 @@
#! FIELDS time iv b.bias b.force2
- 0.000000 1.0485 0.4455 8.9106
- 1.000000 0.9842 0.2741 5.4827
- 2.000000 1.0049 0.3249 6.4986
- 3.000000 0.9474 0.1949 3.8973
- 4.000000 0.9911 0.2907 5.8143
+ 0.000000 1.0002 0.3131 6.2624
+ 1.000000 0.9990 0.3100 6.2005
+ 2.000000 0.9896 0.2870 5.7398
+ 3.000000 1.0004 0.3135 6.2707
+ 4.000000 0.9985 0.3089 6.1773
diff --git a/regtest/crystallization/rt-spherical-integral/forces.reference b/regtest/crystallization/rt-spherical-integral/forces.reference
index 1d0ff64ec..c39663bf3 100644
--- a/regtest/crystallization/rt-spherical-integral/forces.reference
+++ b/regtest/crystallization/rt-spherical-integral/forces.reference
@@ -1,10 +1,10 @@
201
- -0.1864 0.0267 0.3014
-X 0.0003 -0.0003 -0.0003
-X -0.1329 0.1134 -0.1861
-X 0.1329 -0.1134 0.1861
-X 0.0714 0.0810 0.0905
-X -0.0717 -0.0808 -0.0902
+ -0.0027 0.0026 0.0014
+X -0.0000 0.0000 0.0000
+X -0.0018 -0.0021 0.0002
+X 0.0018 0.0021 -0.0002
+X 0.0008 0.0006 0.0009
+X -0.0008 -0.0006 -0.0009
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
@@ -202,12 +202,12 @@ X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
201
- 0.0172 0.1320 -0.1770
-X -0.0001 -0.0000 0.0000
-X 0.0075 -0.0076 -0.0156
-X -0.0310 -0.0846 0.0847
-X -0.0074 0.0077 0.0156
-X 0.0310 0.0846 -0.0847
+ 0.0000 -0.0036 0.0028
+X -0.0000 -0.0000 0.0000
+X 0.0002 -0.0001 -0.0005
+X -0.0003 0.0021 -0.0025
+X -0.0002 0.0001 0.0005
+X 0.0003 -0.0021 0.0025
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
@@ -405,12 +405,12 @@ X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
201
- 0.0128 -0.0660 0.0610
-X 0.0030 0.0044 0.0010
-X -0.0172 0.0225 0.0312
-X 0.0002 -0.0006 -0.0012
-X -0.0015 -0.0029 0.0038
-X 0.0156 -0.0235 -0.0348
+ 0.0125 -0.0616 0.0512
+X 0.0070 0.0067 -0.0013
+X -0.0114 0.0116 0.0292
+X -0.0000 0.0001 0.0002
+X -0.0024 -0.0077 0.0035
+X 0.0069 -0.0107 -0.0316
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
@@ -608,12 +608,12 @@ X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
201
- -0.0155 0.0614 -0.0330
-X 0.0011 -0.0024 -0.0009
-X -0.0060 -0.0023 -0.0019
-X -0.0015 0.0069 0.0035
-X -0.0097 0.1341 0.0065
-X 0.0161 -0.1364 -0.0072
+ -0.0001 -0.0073 0.0092
+X 0.0003 -0.0004 -0.0002
+X 0.0005 -0.0018 -0.0011
+X -0.0006 0.0018 0.0011
+X -0.0013 0.0014 -0.0038
+X 0.0011 -0.0009 0.0040
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
@@ -811,12 +811,12 @@ X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
201
- 0.0124 -0.0030 -0.0094
-X -0.0208 0.0213 0.0286
-X 0.0208 -0.0213 -0.0286
-X -0.0261 0.0137 0.0262
+ 0.0057 -0.0117 0.0061
+X 0.0005 0.0019 0.0006
+X -0.0005 -0.0019 -0.0006
+X 0.0148 0.0097 0.0061
X -0.0000 -0.0000 0.0000
-X 0.0261 -0.0137 -0.0262
+X -0.0148 -0.0097 -0.0061
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
X 0.0000 0.0000 0.0000
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index 3b89af0f2..4b15dee2a 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -111,11 +111,47 @@ void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vec
}
void GridVessel::setupFibonacciGrid( const unsigned& np ) {
- bounds_set=true;
- npoints = np; fib_increment = pi*( 3 - sqrt(5) );
+ bounds_set=true; root5 = sqrt(5);
+ npoints = np; golden = ( 1 + sqrt(5) ) / 2.0; igolden = golden - 1;
+ fib_increment = 2*pi*igolden; log_golden2 = std::log( golden*golden );
fib_offset = 2 / static_cast<double>( npoints );
- Random random; fib_rnd = std::floor( npoints*random.RandU01() );
+ fib_shift = fib_offset/2 - 1;
resize();
+
+ std::vector<double> icoord( dimension ), jcoord( dimension );
+ // Find minimum distance between each pair of points
+ std::vector<unsigned> tindices( dimension ); std::vector<double> mindists( npoints );
+ for(unsigned i=0; i<npoints; ++i) {
+ getFibonacciCoordinates( i, icoord ); mindists[i] = 0;
+ for(unsigned j=0; j<npoints; ++j) {
+ if( i==j ) continue ; // Points are not neighbors to themselves
+ getFibonacciCoordinates( j, jcoord );
+ // Calculate the dot product
+ double dot=0; for(unsigned k=0; k<dimension; ++k) dot += icoord[k]*jcoord[k];
+ if( dot>mindists[i] ) mindists[i]=dot;
+ }
+ }
+ // And now take minimum of dot products
+ double min=mindists[0];
+ for(unsigned i=1; i<npoints; ++i) {
+ if( mindists[i]<min ) min=mindists[i];
+ }
+ double final_cutoff;
+ if( getFibonacciCutoff()<-1 ) final_cutoff=-1;
+ else final_cutoff = cos( acos( getFibonacciCutoff() ) + acos( min ) );
+
+ // And now construct the neighbor list
+ fib_nlist.resize( npoints );
+ for(unsigned i=0; i<npoints; ++i) {
+ getFibonacciCoordinates( i, icoord );
+ for(unsigned j=0; j<npoints; ++j) {
+ if( i==j ) continue ; // Points are not neighbors to themselves
+ getFibonacciCoordinates( j, jcoord );
+ // Calculate the dot product
+ double dot=0; for(unsigned k=0; k<dimension; ++k) dot += icoord[k]*jcoord[k];
+ if( dot>final_cutoff ) { fib_nlist[i].push_back(j); }
+ }
+ }
}
std::string GridVessel::description() {
@@ -170,8 +206,43 @@ void GridVessel::getIndices( const std::vector<double>& point, std::vector<unsig
unsigned GridVessel::getIndex( const std::vector<double>& point ) const {
plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension );
- std::vector<unsigned> indices(dimension); getIndices( point, indices );
- return getIndex( indices );
+ if( gtype==flat ) {
+ std::vector<unsigned> indices(dimension); getIndices( point, indices );
+ return getIndex( indices );
+ } else if( gtype==fibonacci ) {
+ return getFibonacciIndex( point );
+ } else {
+ plumed_error();
+ }
+}
+
+unsigned GridVessel::getFibonacciIndex( const std::vector<double>& p ) const {
+ plumed_dbg_assert( gtype==fibonacci );
+ // Convert input point to coordinates on cylinder
+ unsigned k=2; double phi = atan2( p[2], p[0] ), sinthet2 = 1 - p[1]*p[1];
+ // Calculate power to raise golden ratio
+ if( sinthet2<epsilon ) { k = 2; }
+ else {
+ k = std::floor( std::log( npoints*pi*root5*sinthet2 ) / log_golden2 );
+ if( k<2 ) k = 2;
+ }
+ double Fk = pow( golden, k ) / root5, F0 = std::round(Fk), F1 = std::round(Fk*golden);
+ Matrix<double> B(2,2), invB(2,2); std::vector<double> thisp(3);
+ B(0,0) = 2*pi*((F0+1)*igolden - std::floor((F0+1)*igolden)) - fib_increment;
+ B(0,1) = 2*pi*((F1+1)*igolden - std::floor((F1+1)*igolden)) - fib_increment;
+ B(1,0) = -2*F0/npoints; B(1,1) = -2*F1/npoints; Invert( B, invB );
+ std::vector<double> vv(2), rc(2); vv[0]=-phi; vv[1] = p[1] - fib_shift;
+ mult( invB, vv, rc ); std::vector<int> c(2); c[0]=std::floor(rc[0]); c[1]=std::floor(rc[1]);
+ unsigned outind; double mindist = 10000000.;
+ for(int s=0; s<4; ++s) {
+ double ttt, costheta = B(1,0)*( c[0] + s%2 ) + B(1,1)*( c[1] + s/2 ) + fib_shift;
+ if( costheta>1 ) ttt=1; else if( costheta<-1 ) ttt=-1; else ttt=costheta;
+ costheta = 2*ttt - costheta;
+ unsigned i = std::floor( 0.5*npoints*(1+costheta) ); getFibonacciCoordinates( i, thisp );
+ double dist=0; for(unsigned j=0; j<3; ++j) { double tmp=thisp[j]-p[j]; dist += tmp*tmp; }
+ if( dist<mindist ) { outind = i; mindist = dist; }
+ }
+ return outind;
}
void GridVessel::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
@@ -183,7 +254,6 @@ void GridVessel::convertIndexToIndices( const unsigned& index, const std::vector
if(dimension>=2) { // I think this is wrong
indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
}
-
}
void GridVessel::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
@@ -197,22 +267,27 @@ void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<do
void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
if( gtype==flat ) {
- getIndices( ipoint, tindices );
- for(unsigned i=0; i<dimension; ++i) x[i] = min[i] + dx[i]*tindices[i];
+ getFlatGridCoordinates( ipoint, tindices, x );
} else if( gtype==fibonacci ) {
- x[1] = ((ipoint*fib_offset) - 1) + (fib_offset/2);
- double r = sqrt( 1 -x[1]*x[1] );
- double phi = ((ipoint+fib_rnd)%npoints)*fib_increment;
- x[0] = r*cos(phi); x[2] = r*sin(phi);
- double norm=0;
- for(unsigned j=0; j<3; ++j) norm+=x[j]*x[j];
- norm = sqrt(norm);
- for(unsigned j=0; j<3; ++j) x[j] = x[j] / norm;
+ getFibonacciCoordinates( ipoint, x );
} else {
plumed_error();
}
}
+void GridVessel::getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
+ plumed_dbg_assert( gtype==flat ); getIndices( ipoint, tindices );
+ for(unsigned i=0; i<dimension; ++i) x[i] = min[i] + dx[i]*tindices[i];
+}
+
+void GridVessel::getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
+ plumed_dbg_assert( gtype==fibonacci );
+ x[1] = (ipoint*fib_offset) + fib_shift; double r = sqrt( 1 - x[1]*x[1] );
+ double phi = ipoint*fib_increment; x[0] = r*cos(phi); x[2] = r*sin(phi);
+ double norm=0; for(unsigned j=0; j<3; ++j) norm+=x[j]*x[j];
+ norm = sqrt(norm); for(unsigned j=0; j<3; ++j) x[j] = x[j] / norm;
+}
+
void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const {
plumed_dbg_assert( gtype==flat ); mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
@@ -285,11 +360,20 @@ std::vector<unsigned> GridVessel::getNbin() const {
void GridVessel::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
- plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
-
- std::vector<unsigned> indices( dimension );
- for(unsigned i=0; i<dimension; ++i) indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
- getNeighbors( indices, nneigh, num_neighbors, neighbors );
+ plumed_dbg_assert( bounds_set && nneigh.size()==dimension );
+
+ if( gtype == flat ) {
+ std::vector<unsigned> indices( dimension );
+ for(unsigned i=0; i<dimension; ++i) indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
+ getNeighbors( indices, nneigh, num_neighbors, neighbors );
+ } else if( gtype == fibonacci ) {
+ unsigned find = getFibonacciIndex( pp );
+ num_neighbors = 1 + fib_nlist[find].size();
+ if( neighbors.size()<num_neighbors ) neighbors.resize( num_neighbors );
+ neighbors[0]=find; for(unsigned i=0; i<fib_nlist[find].size(); ++i) neighbors[1+i] = fib_nlist[find][i];
+ } else {
+ plumed_error();
+ }
}
void GridVessel::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
@@ -360,7 +444,7 @@ double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const u
std::vector<unsigned> nindices(dimension);
std::vector<unsigned> indices(dimension); getIndices( x, indices );
std::vector<unsigned> neigh; getSplineNeighbors( getIndex(indices), neigh );
- std::vector<double> xfloor(dimension); getGridPointCoordinates( getIndex(x), xfloor );
+ std::vector<double> xfloor(dimension); getFlatGridCoordinates( getIndex(x), nindices, xfloor );
// loop over neighbors
for(unsigned int ipoint=0; ipoint<neigh.size(); ++ipoint) {
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index 120b9bba4..726ebd9bd 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -41,8 +41,10 @@ private:
/// The number of points in the grid
unsigned npoints;
/// Stuff for fibonacci grids
- int fib_rnd;
- double fib_offset, fib_increment;
+ double root5, golden, igolden, log_golden2;
+/// Fib increment here is equal to 2*pi*(INVERSE GOLDEN RATIO)
+ double fib_offset, fib_increment, fib_shift;
+ std::vector<std::vector<unsigned> > fib_nlist;
/// Units for Gaussian Cube file
double cube_units;
/// This flag is used to check if the user has created a valid input
@@ -81,6 +83,12 @@ protected:
std::vector<bool> active;
/// Convert a point in space the the correspoinding grid point
unsigned getIndex( const std::vector<double>& p ) const ;
+/// Get the index of the closest point on the fibonacci sphere
+ unsigned getFibonacciIndex( const std::vector<double>& p ) const ;
+/// Get the flat grid coordinates
+ void getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const ;
+/// Get the coordinates on the Fibonacci grid
+ void getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const ;
public:
/// keywords
static void registerKeywords( Keywords& keys );
@@ -92,6 +100,8 @@ public:
std::string getType() const ;
/// Set the minimum and maximum of the grid
virtual void setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax, const std::vector<unsigned>& nbins, const std::vector<double>& spacing );
+/// Get the cutoff to use for the Fibonacci spheres
+ virtual double getFibonacciCutoff() const ;
/// Setup the grid if it is a fibonacci grid on the surface of a sphere
void setupFibonacciGrid( const unsigned& np );
/// Get a description of the grid to output to the log
@@ -266,6 +276,11 @@ std::string GridVessel::getType() const {
plumed_error();
}
+inline
+double GridVessel::getFibonacciCutoff() const {
+ return 0.0;
+}
+
}
}
#endif
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index 1ddae05a4..73ec78ab6 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -52,6 +52,10 @@ HistogramOnGrid::HistogramOnGrid( const vesselbase::VesselOptions& da ):
}
}
+double HistogramOnGrid::getFibonacciCutoff() const {
+ return std::log( epsilon / von_misses_norm ) / von_misses_concentration;
+}
+
bool HistogramOnGrid::noDiscreteKernels() const {
return !discrete;
}
@@ -79,10 +83,11 @@ KernelFunctions* HistogramOnGrid::getKernelAndNeighbors( std::vector<double>& po
KernelFunctions* kernel = new KernelFunctions( point, bandwidths, kerneltype, false, 1.0, true );
getNeighbors( kernel->getCenter(), nneigh, num_neigh, neighbors );
return kernel;
+ } else if( getType()=="fibonacci" ) {
+ getNeighbors( point, nneigh, num_neigh, neighbors );
+ return NULL;
} else {
- num_neigh = getNumberOfPoints();
- if( neighbors.size()!=getNumberOfPoints() ) neighbors.resize( getNumberOfPoints() );
- for(unsigned i=0; i<getNumberOfPoints(); ++i) neighbors[i]=i;
+ plumed_error();
}
return NULL;
}
diff --git a/src/gridtools/HistogramOnGrid.h b/src/gridtools/HistogramOnGrid.h
index 96044ffb7..5e4768fef 100644
--- a/src/gridtools/HistogramOnGrid.h
+++ b/src/gridtools/HistogramOnGrid.h
@@ -56,6 +56,7 @@ public:
void addOneKernelEachTimeOnly() { addOneKernelAtATime=true; }
virtual void getFinalForces( const std::vector<double>& buffer, std::vector<double>& finalForces );
bool noDiscreteKernels() const ;
+ double getFibonacciCutoff() const ;
};
inline
--
GitLab