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 : {
46 134 : std::string geom; parse("TYPE",geom);
47 67 : if( geom=="flat" ) gtype=flat;
48 3 : else if( geom=="fibonacci" ) gtype=fibonacci;
49 0 : else plumed_merror( geom + " is invalid geometry type");
50 134 : std::vector<std::string> compnames; parseVector("COMPONENTS",compnames);
51 67 : std::vector<std::string> coordnames; parseVector("COORDINATES",coordnames);
52 67 : if( gtype==flat ) {
53 64 : dimension=coordnames.size();
54 64 : str_min.resize( dimension); str_max.resize( dimension ); stride.resize( dimension );
55 64 : max.resize( dimension ); dx.resize( dimension ); nbin.resize( dimension ); min.resize( dimension );
56 3 : } else if( gtype==fibonacci ) {
57 3 : if( coordnames.size()!=3 ) error("cannot generate fibonacci grid points on surface of sphere if not 3 input coordinates");
58 3 : dimension=3;
59 : }
60 :
61 67 : unsigned n=0; nper=compnames.size()*( 1 + coordnames.size() );
62 67 : arg_names.resize( coordnames.size() + compnames.size()*( 1 + coordnames.size() ) );
63 169 : for(unsigned i=0; i<coordnames.size(); ++i) { arg_names[n] = coordnames[i]; n++; }
64 135 : for(unsigned i=0; i<compnames.size(); ++i) {
65 68 : arg_names[n]=compnames[i]; n++;
66 172 : for(unsigned j=0; j<coordnames.size(); ++j) { arg_names[n] = "d" + compnames[i] + "_" + coordnames[j]; n++; }
67 : }
68 67 : pbc.resize( dimension );
69 67 : std::vector<std::string> spbc( dimension ); parseVector("PBC",spbc);
70 169 : for(unsigned i=0; i<dimension; ++i) {
71 102 : if( spbc[i]=="F" ) pbc[i]=false;
72 31 : else if( spbc[i]=="T" ) pbc[i]=true;
73 0 : else plumed_error();
74 : }
75 134 : }
76 :
77 19 : void GridVessel::setNoDerivatives() {
78 19 : nper = ( nper/(1+dimension) ); noderiv=true;
79 19 : std::vector<std::string> tnames( dimension ), cnames(nper);
80 41 : for(unsigned i=0; i<dimension; ++i) tnames[i]=arg_names[i];
81 38 : unsigned k=dimension; for(unsigned i=0; i<nper; ++i) { cnames[i]=arg_names[k]; k+=(1+dimension); }
82 19 : arg_names.resize( dimension + nper );
83 41 : for(unsigned i=0; i<dimension; ++i) arg_names[i]=tnames[i];
84 38 : for(unsigned i=0; i<nper; ++i) arg_names[dimension+i]=cnames[i];
85 19 : }
86 :
87 98 : void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
88 : const std::vector<unsigned>& binsin, const std::vector<double>& spacing ) {
89 : plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
90 98 : plumed_assert( gtype==flat && (spacing.size()==dimension || binsin.size()==dimension) );
91 :
92 98 : npoints=1; bounds_set=true;
93 245 : for(unsigned i=0; i<dimension; ++i) {
94 147 : str_min[i]=smin[i]; str_max[i]=smax[i];
95 : // GB: I ignore the result here not to break test multicolvar/rt-dens-average
96 : // where these functions were called with str_min[i] and str_max[i] as empty string
97 : // To be checked.
98 147 : Tools::convertNoexcept( str_min[i], min[i] );
99 147 : Tools::convertNoexcept( str_max[i], max[i] );
100 :
101 147 : if( spacing.size()==dimension && binsin.size()==dimension ) {
102 39 : if( spacing[i]==0 ) nbin[i] = binsin[i];
103 : else {
104 37 : double range = max[i] - min[i]; nbin[i] = std::round( range / spacing[i]);
105 : // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
106 37 : if( nbin[i]!=binsin[i] ) plumed_merror("mismatch between input spacing and input number of bins");
107 : }
108 108 : } else if( binsin.size()==dimension ) nbin[i]=binsin[i];
109 0 : else if( spacing.size()==dimension ) nbin[i] = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
110 0 : else plumed_error();
111 147 : dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
112 147 : if( !pbc[i] ) { max[i] +=dx[i]; nbin[i]+=1; }
113 147 : stride[i]=npoints;
114 147 : npoints*=nbin[i];
115 : }
116 98 : resize(); // Always resize after setting new bounds as grid size may have have changed
117 98 : }
118 :
119 3 : void GridVessel::setupFibonacciGrid( const unsigned& np ) {
120 3 : bounds_set=true; root5 = std::sqrt(5);
121 3 : npoints = np; golden = ( 1 + std::sqrt(5) ) / 2.0; igolden = golden - 1;
122 3 : fib_increment = 2*pi*igolden; log_golden2 = std::log( golden*golden );
123 3 : fib_offset = 2 / static_cast<double>( npoints );
124 3 : fib_shift = fib_offset/2 - 1;
125 3 : resize();
126 :
127 3 : std::vector<double> icoord( dimension ), jcoord( dimension );
128 : // Find minimum distance between each pair of points
129 3 : std::vector<double> mindists( npoints );
130 347 : for(unsigned i=0; i<npoints; ++i) {
131 344 : getFibonacciCoordinates( i, icoord ); mindists[i] = 0;
132 41080 : for(unsigned j=0; j<npoints; ++j) {
133 40736 : if( i==j ) continue ; // Points are not neighbors to themselves
134 40392 : getFibonacciCoordinates( j, jcoord );
135 : // Calculate the dot product
136 161568 : double dot=0; for(unsigned k=0; k<dimension; ++k) dot += icoord[k]*jcoord[k];
137 40392 : if( dot>mindists[i] ) mindists[i]=dot;
138 : }
139 : }
140 : // And now take minimum of dot products
141 3 : double min=mindists[0];
142 344 : for(unsigned i=1; i<npoints; ++i) {
143 341 : if( mindists[i]<min ) min=mindists[i];
144 : }
145 : double final_cutoff;
146 3 : if( getFibonacciCutoff()<-1 ) final_cutoff=-1;
147 2 : else final_cutoff = std::cos( std::acos( getFibonacciCutoff() ) + std::acos( min ) );
148 :
149 : // And now construct the neighbor list
150 3 : fib_nlist.resize( npoints );
151 347 : for(unsigned i=0; i<npoints; ++i) {
152 344 : getFibonacciCoordinates( i, icoord );
153 41080 : for(unsigned j=0; j<npoints; ++j) {
154 40736 : if( i==j ) continue ; // Points are not neighbors to themselves
155 40392 : getFibonacciCoordinates( j, jcoord );
156 : // Calculate the dot product
157 161568 : double dot=0; for(unsigned k=0; k<dimension; ++k) dot += icoord[k]*jcoord[k];
158 40392 : if( dot>final_cutoff ) { fib_nlist[i].push_back(j); }
159 : }
160 : }
161 3 : }
162 :
163 65 : std::string GridVessel::description() {
164 65 : if( !bounds_set ) return "";
165 :
166 : std::string des;
167 51 : if( gtype==flat ) {
168 : des="grid of "; std::string num;
169 68 : for(unsigned i=0; i<dimension-1; ++i) {
170 19 : Tools::convert( nbin[i], num );
171 38 : des += num + " X ";
172 : }
173 49 : Tools::convert( nbin[dimension-1], num );
174 49 : des += num + " equally spaced points between (";
175 68 : for(unsigned i=0; i<dimension-1; ++i) des += str_min[i] + ",";
176 49 : Tools::convert( nbin[dimension-1], num );
177 49 : des += str_min[dimension-1] + ") and (";
178 68 : for(unsigned i=0; i<dimension-1; ++i) des += str_max[i] + ",";
179 98 : des += str_max[dimension-1] + ")";
180 2 : } else if( gtype==fibonacci ) {
181 2 : std::string num; Tools::convert( npoints, num );
182 4 : des += "fibonacci grid of " + num + " points on spherical surface";
183 : }
184 : return des;
185 : }
186 :
187 213 : void GridVessel::resize() {
188 213 : plumed_massert( nper>0, "Number of datapoints at each grid point has not been set");
189 213 : if( getAction() ) resizeBuffer( getNumberOfBufferPoints()*nper + 1 + 2*getAction()->getNumberOfDerivatives() );
190 213 : setDataSize( npoints*nper ); forces.resize( npoints );
191 213 : if( active.size()!=npoints) active.resize( npoints, true );
192 213 : }
193 :
194 24333356 : unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
195 : plumed_dbg_assert( gtype==flat && bounds_set && indices.size()==dimension );
196 : // indices are flattended using a column-major order
197 24333356 : unsigned index=indices[dimension-1];
198 69887354 : for(unsigned i=dimension-1; i>0; --i) {
199 45553998 : index=index*nbin[i-1]+indices[i-1];
200 : }
201 24333356 : return index;
202 : }
203 :
204 1060385 : void GridVessel::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
205 : plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension && indices.size()==dimension );
206 3196118 : for(unsigned i=0; i<dimension; ++i) {
207 2135733 : indices[i]=std::floor( (point[i] - min[i])/dx[i] );
208 2135733 : if( pbc[i] ) indices[i]=indices[i]%nbin[i];
209 2090221 : else if( indices[i]>nbin[i] ) {
210 0 : std::string pp, num; Tools::convert( point[0], pp );
211 0 : for(unsigned j=1; j<point.size(); ++j) { Tools::convert( point[j], num ); pp += ", " + num; }
212 0 : plumed_merror("point (" + pp + ") is outside grid range");
213 : }
214 : }
215 1060385 : }
216 :
217 530954 : unsigned GridVessel::getIndex( const std::vector<double>& point ) const {
218 : plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension );
219 530954 : if( gtype==flat ) {
220 530954 : std::vector<unsigned> indices(dimension); getIndices( point, indices );
221 530954 : return getIndex( indices );
222 0 : } else if( gtype==fibonacci ) {
223 0 : return getFibonacciIndex( point );
224 : } else {
225 0 : plumed_error();
226 : }
227 : }
228 :
229 57 : unsigned GridVessel::getFibonacciIndex( const std::vector<double>& p ) const {
230 : plumed_dbg_assert( gtype==fibonacci );
231 : // Convert input point to coordinates on cylinder
232 57 : int k=2; double phi = std::atan2( p[2], p[0] ), sinthet2 = 1 - p[1]*p[1];
233 : // Calculate power to raise golden ratio
234 57 : if( sinthet2<epsilon ) { k = 2; }
235 : else {
236 57 : k = std::floor( std::log( npoints*pi*root5*sinthet2 ) / log_golden2 );
237 : if( k<2 ) k = 2;
238 : }
239 57 : double Fk = pow( golden, k ) / root5, F0 = std::round(Fk), F1 = std::round(Fk*golden);
240 57 : Matrix<double> B(2,2), invB(2,2); std::vector<double> thisp(3);
241 57 : B(0,0) = 2*pi*((F0+1)*igolden - std::floor((F0+1)*igolden)) - fib_increment;
242 57 : B(0,1) = 2*pi*((F1+1)*igolden - std::floor((F1+1)*igolden)) - fib_increment;
243 57 : B(1,0) = -2*F0/npoints; B(1,1) = -2*F1/npoints; Invert( B, invB );
244 57 : std::vector<double> vv(2), rc(2); vv[0]=-phi; vv[1] = p[1] - fib_shift;
245 57 : mult( invB, vv, rc ); std::vector<int> c(2); c[0]=std::floor(rc[0]); c[1]=std::floor(rc[1]);
246 : unsigned outind=0; double mindist = 10000000.;
247 285 : for(int s=0; s<4; ++s) {
248 228 : double ttt, costheta = B(1,0)*( c[0] + s%2 ) + B(1,1)*( c[1] + s/2 ) + fib_shift;
249 228 : if( costheta>1 ) ttt=1; else if( costheta<-1 ) ttt=-1; else ttt=costheta;
250 228 : costheta = 2*ttt - costheta;
251 228 : unsigned i = std::floor( 0.5*npoints*(1+costheta) ); getFibonacciCoordinates( i, thisp );
252 912 : double dist=0; for(unsigned j=0; j<3; ++j) { double tmp=thisp[j]-p[j]; dist += tmp*tmp; }
253 228 : if( dist<mindist ) { outind = i; mindist = dist; }
254 : }
255 57 : return outind;
256 : }
257 :
258 112427797 : void GridVessel::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
259 112427797 : plumed_dbg_assert( gtype==flat ); unsigned kk=index; indices[0]=index%nnbin[0];
260 220890721 : for(unsigned i=1; i<dimension-1; ++i) {
261 108462924 : kk=(kk-indices[i-1])/nnbin[i-1];
262 108462924 : indices[i]=kk%nnbin[i];
263 : }
264 112427797 : if(dimension>=2) { // I think this is wrong
265 112411290 : indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
266 : }
267 112427797 : }
268 :
269 15878999 : void GridVessel::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
270 15878999 : plumed_dbg_assert( gtype==flat ); convertIndexToIndices( index, nbin, indices );
271 15878999 : }
272 :
273 665380 : void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
274 665380 : std::vector<unsigned> tindices( dimension ); getGridPointCoordinates( ipoint, tindices, x );
275 665380 : }
276 :
277 12514073 : void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
278 : plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
279 12514073 : if( gtype==flat ) {
280 12508956 : getFlatGridCoordinates( ipoint, tindices, x );
281 5117 : } else if( gtype==fibonacci ) {
282 5117 : getFibonacciCoordinates( ipoint, x );
283 : } else {
284 0 : plumed_error();
285 : }
286 12514073 : }
287 :
288 13037799 : void GridVessel::getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
289 13037799 : plumed_dbg_assert( gtype==flat ); getIndices( ipoint, tindices );
290 50950437 : for(unsigned i=0; i<dimension; ++i) x[i] = min[i] + dx[i]*tindices[i];
291 13037799 : }
292 :
293 86817 : void GridVessel::getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
294 : plumed_dbg_assert( gtype==fibonacci );
295 86817 : x[1] = (ipoint*fib_offset) + fib_shift; double r = std::sqrt( 1 - x[1]*x[1] );
296 86817 : double phi = ipoint*fib_increment; x[0] = r*std::cos(phi); x[2] = r*std::sin(phi);
297 347268 : double norm=0; for(unsigned j=0; j<3; ++j) norm+=x[j]*x[j];
298 347268 : norm = std::sqrt(norm); for(unsigned j=0; j<3; ++j) x[j] = x[j] / norm;
299 86817 : }
300 :
301 528843 : void GridVessel::getSplineNeighbors( const unsigned& mybox, unsigned& nneighbors, std::vector<unsigned>& mysneigh ) const {
302 528843 : plumed_dbg_assert( gtype==flat ); unsigned nneigh=unsigned(pow(2.0,int(dimension)));
303 528843 : if( mysneigh.size()!=nneigh ) mysneigh.resize(nneigh);
304 :
305 528843 : unsigned inind; nneighbors = 0;
306 528843 : std::vector<unsigned> tmp_indices( dimension );
307 528843 : std::vector<unsigned> my_indices( dimension );
308 528843 : getIndices( mybox, my_indices );
309 2677587 : for(unsigned i=0; i<nneigh; ++i) {
310 : unsigned tmp=i; inind=0;
311 6513376 : for(unsigned j=0; j<dimension; ++j) {
312 4364632 : unsigned i0=tmp%2+my_indices[j]; tmp/=2;
313 4364632 : if(!pbc[j] && i0==nbin[j]) continue;
314 4323832 : if( pbc[j] && i0==nbin[j]) i0=0;
315 4323832 : tmp_indices[inind++]=i0;
316 : }
317 2148744 : if(inind==dimension ) {
318 2108144 : unsigned findex=getIndex( tmp_indices );
319 2108144 : mysneigh[nneighbors++]=findex;
320 2108144 : plumed_massert( active[findex], "inactive grid point required for splines");
321 : }
322 : }
323 528843 : }
324 :
325 6552681 : double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
326 13105362 : plumed_assert( bounds_set && ipoint<npoints && jelement<nper && active[ipoint] );
327 6552681 : return getDataElement( nper*ipoint + jelement );
328 : }
329 :
330 154208 : void GridVessel::setGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ) {
331 : plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
332 154208 : setDataElement( nper*ipoint + jelement, value );
333 154208 : }
334 :
335 24200 : void GridVessel::setValueAndDerivatives( const unsigned& ipoint, const unsigned& jelement, const double& value, const std::vector<double>& der ) {
336 : plumed_dbg_assert( !noderiv && jelement<getNumberOfComponents() && der.size()==nbin.size() );
337 72600 : setGridElement( ipoint, jelement, value ); for(unsigned i=0; i<der.size(); ++i) setGridElement( ipoint, jelement+1+i, der[i] );
338 24200 : }
339 :
340 2605 : void GridVessel::addToGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ) {
341 : plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
342 2605 : addDataElement( nper*ipoint + jelement, value );
343 2605 : }
344 :
345 15087 : void GridVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
346 : plumed_dbg_assert( myvals.getNumberOfValues()==(nper+1) );
347 65138 : for(unsigned i=0; i<nper; ++i) buffer[bufstart + nper*current + i] += myvals.get(i+1);
348 15087 : }
349 :
350 118 : void GridVessel::finish( const std::vector<double>& buffer ) {
351 118 : if( wasforced ) getFinalForces( buffer, finalForces );
352 98 : else AveragingVessel::finish( buffer );
353 118 : }
354 :
355 0 : double GridVessel::getGridElement( const std::vector<unsigned>& indices, const unsigned& jelement ) const {
356 0 : return getGridElement( getIndex( indices ), jelement );
357 : }
358 :
359 0 : void GridVessel::setGridElement( const std::vector<unsigned>& indices, const unsigned& jelement, const double& value ) {
360 0 : setGridElement( getIndex( indices ), jelement, value );
361 0 : }
362 :
363 202520 : std::vector<std::string> GridVessel::getMin() const {
364 202520 : plumed_dbg_assert( gtype==flat ); return str_min;
365 : }
366 :
367 202553 : std::vector<std::string> GridVessel::getMax() const {
368 202553 : plumed_dbg_assert( gtype==flat ); return str_max;
369 : }
370 :
371 203150 : std::vector<unsigned> GridVessel::getNbin() const {
372 : plumed_dbg_assert( gtype==flat && bounds_set );
373 203150 : std::vector<unsigned> ngrid( dimension );
374 605377 : for(unsigned i=0; i<dimension; ++i) {
375 402227 : if( !pbc[i] ) ngrid[i]=nbin[i] - 1;
376 34200 : else ngrid[i]=nbin[i];
377 : }
378 203150 : return ngrid;
379 : }
380 :
381 21514 : void GridVessel::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
382 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
383 : plumed_dbg_assert( bounds_set );
384 :
385 21514 : if( gtype == flat ) {
386 : plumed_dbg_assert( nneigh.size()==dimension );
387 21457 : std::vector<unsigned> indices( dimension );
388 85340 : for(unsigned i=0; i<dimension; ++i) indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
389 21457 : getNeighbors( indices, nneigh, num_neighbors, neighbors );
390 57 : } else if( gtype == fibonacci ) {
391 57 : unsigned find = getFibonacciIndex( pp );
392 57 : num_neighbors = 1 + fib_nlist[find].size();
393 57 : if( neighbors.size()<num_neighbors ) neighbors.resize( num_neighbors );
394 4773 : neighbors[0]=find; for(unsigned i=0; i<fib_nlist[find].size(); ++i) neighbors[1+i] = fib_nlist[find][i];
395 : } else {
396 0 : plumed_error();
397 : }
398 21514 : }
399 :
400 40878 : void GridVessel::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
401 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
402 : plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
403 :
404 40878 : unsigned num_neigh=1; std::vector<unsigned> small_bin( dimension );
405 163024 : for(unsigned i=0; i<dimension; ++i) {
406 122146 : small_bin[i]=(2*nneigh[i]+1);
407 122146 : num_neigh *=small_bin[i];
408 : }
409 40878 : if( neighbors.size()!=num_neigh ) neighbors.resize( num_neigh );
410 :
411 40878 : num_neighbors=0;
412 40878 : std::vector<unsigned> s_indices(dimension), t_indices(dimension);
413 96589676 : for(unsigned index=0; index<num_neigh; ++index) {
414 : bool found=true;
415 96548798 : convertIndexToIndices( index, small_bin, s_indices );
416 386166232 : for(unsigned i=0; i<dimension; ++i) {
417 289617434 : int i0=s_indices[i]-nneigh[i]+indices[i];
418 289617434 : if(!pbc[i] && i0<0) found=false;
419 289617434 : if(!pbc[i] && i0>=nbin[i]) found=false;
420 289617434 : if( pbc[i] && i0<0) i0=nbin[i]-(-i0)%nbin[i];
421 289617434 : if( pbc[i] && i0>=nbin[i]) i0%=nbin[i];
422 289617434 : t_indices[i]=static_cast<unsigned>(i0);
423 : }
424 96548798 : if( found ) {
425 21093803 : neighbors[num_neighbors]=getIndex( t_indices );
426 21093803 : num_neighbors++;
427 : }
428 : }
429 40878 : }
430 :
431 11 : void GridVessel::setCubeUnits( const double& units ) {
432 11 : plumed_dbg_assert( gtype==flat ); cube_units=units;
433 11 : }
434 :
435 8 : double GridVessel::getCubeUnits() const {
436 8 : plumed_dbg_assert( gtype==flat ); return cube_units;
437 : }
438 :
439 17 : std::string GridVessel::getInputString() const {
440 17 : std::string mstring="COORDINATES="+arg_names[0];
441 23 : for(unsigned i=1; i<dimension; ++i) mstring+="," + arg_names[i];
442 17 : if( gtype==flat ) {
443 : mstring += " TYPE=flat PBC=";
444 17 : if( pbc[0] ) mstring +="T";
445 : else mstring +="F";
446 23 : for(unsigned i=1; i<dimension; ++i) {
447 6 : if( pbc[i] ) mstring +=",T";
448 : else mstring +=",F";
449 : }
450 0 : } else if( gtype==fibonacci ) {
451 : mstring += " TYPE=fibonacci";
452 : }
453 17 : return mstring;
454 : }
455 :
456 528843 : double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const {
457 : plumed_dbg_assert( gtype==flat && der.size()==dimension && !noderiv && ind<getNumberOfComponents() );
458 :
459 528843 : double X,X2,X3,value=0; der.assign(der.size(),0.0);
460 528843 : std::vector<double> fd(dimension);
461 528843 : std::vector<double> C(dimension);
462 528843 : std::vector<double> D(dimension);
463 528843 : std::vector<double> dder(dimension);
464 :
465 528843 : std::vector<unsigned> nindices(dimension); unsigned n_neigh;
466 528843 : std::vector<unsigned> indices(dimension); getIndices( x, indices );
467 528843 : std::vector<unsigned> neigh; getSplineNeighbors( getIndex(indices), n_neigh, neigh );
468 528843 : std::vector<double> xfloor(dimension); getFlatGridCoordinates( getIndex(x), nindices, xfloor );
469 :
470 : // loop over neighbors
471 2636987 : for(unsigned int ipoint=0; ipoint<n_neigh; ++ipoint) {
472 2108144 : double grid=getGridElement(neigh[ipoint], ind*(1+dimension) );
473 6391576 : for(unsigned j=0; j<dimension; ++j) dder[j] = getGridElement( neigh[ipoint], ind*(1+dimension) + 1 + j );
474 :
475 2108144 : getIndices( neigh[ipoint], nindices );
476 : double ff=1.0;
477 :
478 6391576 : for(unsigned j=0; j<dimension; ++j) {
479 : int x0=1;
480 4283432 : if(nindices[j]==indices[j]) x0=0;
481 4283432 : double ddx=dx[j];
482 4283432 : X=std::fabs((x[j]-xfloor[j])/ddx-(double)x0);
483 4283432 : X2=X*X;
484 4283432 : X3=X2*X;
485 : double yy;
486 4283432 : if(std::fabs(grid)<0.0000001) yy=0.0;
487 4269585 : else yy=-dder[j]/grid;
488 6445348 : C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*ddx;
489 4283432 : D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*ddx;
490 4283432 : D[j]*=(x0?-1.0:1.0)/ddx;
491 4283432 : ff*=C[j];
492 : }
493 6391576 : for(unsigned j=0; j<dimension; ++j) {
494 4283432 : fd[j]=D[j];
495 13052528 : for(unsigned i=0; i<dimension; ++i) if(i!=j) fd[j]*=C[i];
496 : }
497 2108144 : value+=grid*ff;
498 6391576 : for(unsigned j=0; j<dimension; ++j) der[j]+=grid*fd[j];
499 : }
500 528843 : return value;
501 : }
502 :
503 3 : void GridVessel::activateThesePoints( const std::vector<bool>& to_activate ) {
504 : plumed_dbg_assert( to_activate.size()==npoints );
505 29991 : for(unsigned i=0; i<npoints; ++i) active[i]=to_activate[i];
506 3 : }
507 :
508 20 : void GridVessel::setForce( const std::vector<double>& inforces ) {
509 : plumed_dbg_assert( inforces.size()==npoints );
510 5725 : wasforced=true; for(unsigned i=0; i<npoints; ++i) forces[i]=inforces[i];
511 20 : }
512 :
513 11870273 : bool GridVessel::wasForced() const {
514 11870273 : return wasforced;
515 : }
516 :
517 20 : bool GridVessel::applyForce( std::vector<double>& fforces ) {
518 : plumed_dbg_assert( fforces.size()==finalForces.size() );
519 20 : if( !wasforced ) return false;
520 3815 : for(unsigned i=0; i<finalForces.size(); ++i) fforces[i]=finalForces[i];
521 20 : wasforced=false; return true;
522 : }
523 :
524 : }
525 : }
526 :
|
