Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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 : #ifndef __PLUMED_gridtools_GridVessel_h
23 : #define __PLUMED_gridtools_GridVessel_h
24 :
25 : #include <string>
26 : #include <cstring>
27 : #include <vector>
28 : #include "vesselbase/AveragingVessel.h"
29 :
30 : namespace PLMD {
31 : namespace gridtools {
32 :
33 199 : class GridVessel : public vesselbase::AveragingVessel {
34 : friend class ActionWithInputGrid;
35 : friend class DumpGrid;
36 : private:
37 : /// The way that grid points are constructed
38 : enum {flat,fibonacci} gtype;
39 : /// Have the minimum and maximum for the grid been set
40 : bool bounds_set;
41 : /// The number of points in the grid
42 : unsigned npoints;
43 : /// Stuff for fibonacci grids
44 : double root5, golden, igolden, log_golden2;
45 : /// Fib increment here is equal to 2*pi*(INVERSE GOLDEN RATIO)
46 : double fib_offset, fib_increment, fib_shift;
47 : std::vector<std::vector<unsigned> > fib_nlist;
48 : /// Units for Gaussian Cube file
49 : double cube_units;
50 : /// This flag is used to check if the user has created a valid input
51 : bool foundprint;
52 : /// The minimum and maximum of the grid stored as doubles
53 : std::vector<double> min, max;
54 : /// The numerical distance between adjacent grid points
55 : std::vector<unsigned> stride;
56 : /// The number of bins in each grid direction
57 : std::vector<unsigned> nbin;
58 : /// The grid point that was requested last by getGridPointCoordinates
59 : unsigned currentGridPoint;
60 : /// The forces that will be output at the end of the calculation
61 : std::vector<double> finalForces;
62 : protected:
63 : /// Is forced
64 : bool wasforced;
65 : /// Forces acting on grid elements
66 : std::vector<double> forces;
67 : /// Do we have derivatives
68 : bool noderiv;
69 : /// The names of the various columns in the grid file
70 : std::vector<std::string> arg_names;
71 : /// The number of pieces of information we are storing for each
72 : /// point in the grid
73 : unsigned nper;
74 : /// Is this direction periodic
75 : std::vector<bool> pbc;
76 : /// The minimum and maximum in the grid stored as strings
77 : std::vector<std::string> str_min, str_max;
78 : /// The spacing between grid points
79 : std::vector<double> dx;
80 : /// The dimensionality of the grid
81 : unsigned dimension;
82 : /// Which grid points are we actively accumulating
83 : std::vector<bool> active;
84 : /// Convert a point in space the the correspoinding grid point
85 : unsigned getIndex( const std::vector<double>& p ) const ;
86 : /// Get the index of the closest point on the fibonacci sphere
87 : unsigned getFibonacciIndex( const std::vector<double>& p ) const ;
88 : /// Get the flat grid coordinates
89 : void getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const ;
90 : /// Get the coordinates on the Fibonacci grid
91 : void getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const ;
92 : public:
93 : /// keywords
94 : static void registerKeywords( Keywords& keys );
95 : /// Constructor
96 : explicit GridVessel( const vesselbase::VesselOptions& );
97 : /// Remove the derivatives
98 : void setNoDerivatives();
99 : /// Get the type of grid we are using
100 : std::string getType() const ;
101 : /// Set the minimum and maximum of the grid
102 : 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 );
103 : /// Get the cutoff to use for the Fibonacci spheres
104 : virtual double getFibonacciCutoff() const ;
105 : /// Setup the grid if it is a fibonacci grid on the surface of a sphere
106 : void setupFibonacciGrid( const unsigned& np );
107 : /// Get a description of the grid to output to the log
108 : std::string description();
109 : /// Convert an index into indices
110 : void convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const ;
111 : /// Flatten the grid and get the grid index for a point
112 : unsigned getIndex( const std::vector<unsigned>& indices ) const ;
113 : /// Get the indices fof a point
114 : void getIndices( const unsigned& index, std::vector<unsigned>& indices ) const ;
115 : /// Get the indices of a particular point
116 : void getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const ;
117 : /// Operations on one of the elements of grid point i
118 : void setGridElement( const unsigned&, const unsigned&, const double& );
119 : /// Add data to an element of the grid
120 : void addToGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value );
121 : /// Operations on one of the elements of grid point specified by vector
122 : double getGridElement( const std::vector<unsigned>&, const unsigned& ) const ;
123 : void setGridElement( const std::vector<unsigned>&, const unsigned&, const double& );
124 : /// Set the size of the buffer equal to nper*npoints
125 : virtual void resize();
126 : /// Get the number of points in the grid
127 : unsigned getNumberOfPoints() const;
128 : /// Get the coordinates for a point in the grid
129 : void getGridPointCoordinates( const unsigned&, std::vector<double>& ) const ;
130 : void getGridPointCoordinates( const unsigned&, std::vector<unsigned>&, std::vector<double>& ) const ;
131 : /// Get the dimensionality of the function
132 : unsigned getDimension() const ;
133 : /// Get the number of components in the vector stored on each grid point
134 : virtual unsigned getNumberOfComponents() const ;
135 : /// Is the grid periodic in the ith direction
136 : bool isPeriodic( const unsigned& i ) const ;
137 : /// Get the number of quantities we have stored at each grid point
138 : unsigned getNumberOfQuantities() const ;
139 : /// Get the number of grid points for each dimension
140 : std::vector<unsigned> getNbin() const ;
141 : /// Get the name of the ith component
142 : std::string getComponentName( const unsigned& i ) const ;
143 : /// Get the vector containing the minimum value of the grid in each dimension
144 : std::vector<std::string> getMin() const ;
145 : /// Get the vector containing the maximum value of the grid in each dimension
146 : std::vector<std::string> getMax() const ;
147 : /// Get the number of points needed in the buffer
148 : virtual unsigned getNumberOfBufferPoints() const ;
149 : /// Get the stride (the distance between the grid points of an index)
150 : const std::vector<unsigned>& getStride() const ;
151 : /// Return the volume of one of the grid cells
152 : double getCellVolume() const ;
153 : /// Get the value of the ith grid element
154 : virtual double getGridElement( const unsigned&, const unsigned& ) const ;
155 : /// Get the set of points neighouring a particular location in space
156 : void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
157 : unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
158 : /// Get the neighbors for a set of indices of a point
159 : void getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
160 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const ;
161 : /// Get the points neighboring a particular spline point
162 : void getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const ;
163 : /// Get the spacing between grid points
164 : const std::vector<double>& getGridSpacing() const ;
165 : /// Get the extent of the grid in one of the axis
166 : double getGridExtent( const unsigned& i ) const ;
167 : /// Copy data from the action into the grid
168 : virtual void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
169 : /// Finish the calculation
170 : virtual void finish( const std::vector<double>& buffer );
171 : /// This ensures that Gaussian cube fies are in correct units
172 : void setCubeUnits( const double& units );
173 : /// This ensures that Gaussian cube files are in correct units
174 : double getCubeUnits() const ;
175 : /// Return a string containing the input to the grid so we can clone it
176 : std::string getInputString() const ;
177 : /// Does this have derivatives
178 : bool noDerivatives() const ;
179 : /// Get the value and derivatives at a particular location using spline interpolation
180 : double getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const ;
181 : /// Deactivate all the grid points
182 : void activateThesePoints( const std::vector<bool>& to_activate );
183 : /// Is this point active
184 : bool inactive( const unsigned& ip ) const ;
185 : /// This retrieves the final force
186 0 : virtual void getFinalForces( const std::vector<double>& buffer, std::vector<double>& finalForces ) { plumed_error(); }
187 : /// Apply the forces
188 : void setForce( const std::vector<double>& inforces );
189 : /// Was a force added to the grid
190 : bool wasForced() const ;
191 : /// And retrieve the forces
192 : bool applyForce( std::vector<double>& fforces );
193 : };
194 :
195 : inline
196 : unsigned GridVessel::getNumberOfQuantities() const {
197 357499 : return nper;
198 : }
199 :
200 : inline
201 : unsigned GridVessel::getNumberOfPoints() const {
202 195540 : return npoints;
203 : }
204 :
205 : inline
206 599 : const std::vector<double>& GridVessel::getGridSpacing() const {
207 599 : if( gtype==flat ) return dx;
208 0 : plumed_merror("dont understand what spacing means for spherical grids");
209 : return dx;
210 : }
211 :
212 : inline
213 : double GridVessel::getCellVolume() const {
214 4 : if( gtype==flat ) {
215 13 : double myvol=1.0; for(unsigned i=0; i<dimension; ++i) myvol *= dx[i];
216 : return myvol;
217 : } else {
218 1 : return 4*pi / static_cast<double>( getNumberOfPoints() );
219 : }
220 : }
221 :
222 : inline
223 : unsigned GridVessel::getDimension() const {
224 819175 : return dimension;
225 : }
226 :
227 : inline
228 : bool GridVessel::isPeriodic( const unsigned& i ) const {
229 : plumed_dbg_assert( gtype==flat );
230 40787 : return pbc[i];
231 : }
232 :
233 : inline
234 : std::string GridVessel::getComponentName( const unsigned& i ) const {
235 64 : return arg_names[i];
236 : }
237 :
238 : inline
239 30 : unsigned GridVessel::getNumberOfComponents() const {
240 30 : if( noderiv ) return nper;
241 27 : return nper / ( dimension + 1 );
242 : }
243 :
244 : inline
245 : double GridVessel::getGridExtent( const unsigned& i ) const {
246 : plumed_dbg_assert( gtype==flat );
247 92 : return max[i] - min[i];
248 : }
249 :
250 : inline
251 : bool GridVessel::noDerivatives() const {
252 11365 : return noderiv;
253 : }
254 :
255 : inline
256 : bool GridVessel::inactive( const unsigned& ip ) const {
257 : plumed_dbg_assert( ip<npoints );
258 21106479 : return !active[ip];
259 : }
260 :
261 : inline
262 : const std::vector<unsigned>& GridVessel::getStride() const {
263 : plumed_dbg_assert( gtype==flat );
264 : return stride;
265 : }
266 :
267 : inline
268 44 : unsigned GridVessel::getNumberOfBufferPoints() const {
269 277 : return npoints;
270 : }
271 :
272 : inline
273 21958 : std::string GridVessel::getType() const {
274 21958 : if( gtype==flat ) return "flat";
275 422 : else if( gtype==fibonacci ) return "fibonacci";
276 0 : plumed_error();
277 : }
278 :
279 : inline
280 2 : double GridVessel::getFibonacciCutoff() const {
281 2 : return 0.0;
282 : }
283 :
284 : }
285 : }
286 : #endif
|
