Skip to content
Snippets Groups Projects
Commit fd27f122 authored by Gareth Tribello's avatar Gareth Tribello
Browse files

Added routines for truncating kernel evaluation on spherical fibonacci grids

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