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             : #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             : 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() override;
     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 values and derivatives of a particular element
     125             :   void setValueAndDerivatives( const unsigned&, const unsigned&, const double&, const std::vector<double>& );
     126             : /// Set the size of the buffer equal to nper*npoints
     127             :   void resize() override;
     128             : /// Get the number of points in the grid
     129             :   unsigned getNumberOfPoints() const;
     130             : /// Get the coordinates for a point in the grid
     131             :   void getGridPointCoordinates( const unsigned&, std::vector<double>& ) const ;
     132             :   void getGridPointCoordinates( const unsigned&, std::vector<unsigned>&, std::vector<double>& ) const ;
     133             : /// Get the dimensionality of the function
     134             :   unsigned getDimension() const ;
     135             : /// Get the number of components in the vector stored on each grid point
     136             :   virtual unsigned getNumberOfComponents() const ;
     137             : /// Is the grid periodic in the ith direction
     138             :   bool isPeriodic( const unsigned& i ) const ;
     139             : /// Get the number of quantities we have stored at each grid point
     140             :   unsigned getNumberOfQuantities() const ;
     141             : /// Get the number of grid points for each dimension
     142             :   std::vector<unsigned> getNbin() const ;
     143             : /// Get the name of the ith component
     144             :   std::string getComponentName( const unsigned& i ) const ;
     145             : /// Get the vector containing the minimum value of the grid in each dimension
     146             :   std::vector<std::string> getMin() const ;
     147             : /// Get the vector containing the maximum value of the grid in each dimension
     148             :   std::vector<std::string> getMax() const ;
     149             : /// Get the number of points needed in the buffer
     150             :   virtual unsigned getNumberOfBufferPoints() const ;
     151             : /// Get the stride (the distance between the grid points of an index)
     152             :   const std::vector<unsigned>& getStride() const ;
     153             : /// Return the volume of one of the grid cells
     154             :   double getCellVolume() const ;
     155             : /// Get the value of the ith grid element
     156             :   virtual double getGridElement( const unsigned&, const unsigned& ) const ;
     157             : /// Get the set of points neighouring a particular location in space
     158             :   void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
     159             :                      unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
     160             : /// Get the neighbors for a set of indices of a point
     161             :   void getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
     162             :                      unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const ;
     163             : /// Get the points neighboring a particular spline point
     164             :   void getSplineNeighbors( const unsigned& mybox, unsigned& nneighbors, std::vector<unsigned>& mysneigh ) const ;
     165             : /// Get the spacing between grid points
     166             :   const std::vector<double>& getGridSpacing() const ;
     167             : /// Get the extent of the grid in one of the axis
     168             :   double getGridExtent( const unsigned& i ) const ;
     169             : /// Copy data from the action into the grid
     170             :   void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const override;
     171             : /// Finish the calculation
     172             :   void finish( const std::vector<double>& buffer ) override;
     173             : /// This ensures that Gaussian cube fies are in correct units
     174             :   void setCubeUnits( const double& units );
     175             : /// This ensures that Gaussian cube files are in correct units
     176             :   double getCubeUnits() const ;
     177             : /// Return a string containing the input to the grid so we can clone it
     178             :   std::string getInputString() const ;
     179             : /// Does this have derivatives
     180             :   bool noDerivatives() const ;
     181             : /// Get the value and derivatives at a particular location using spline interpolation
     182             :   double getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const ;
     183             : /// Deactivate all the grid points
     184             :   void activateThesePoints( const std::vector<bool>& to_activate );
     185             : /// Is this point active
     186             :   bool inactive( const unsigned& ip ) const ;
     187             : /// This retrieves the final force
     188           0 :   virtual void getFinalForces( const std::vector<double>& buffer, std::vector<double>& finalForces ) { plumed_error(); }
     189             : /// Apply the forces
     190             :   void setForce( const std::vector<double>& inforces );
     191             : /// Was a force added to the grid
     192             :   bool wasForced() const ;
     193             : /// And retrieve the forces
     194             :   bool applyForce( std::vector<double>& fforces ) override;
     195             : };
     196             : 
     197             : inline
     198             : unsigned GridVessel::getNumberOfQuantities() const {
     199      440217 :   return nper;
     200             : }
     201             : 
     202             : inline
     203             : unsigned GridVessel::getNumberOfPoints() const {
     204      873010 :   return npoints;
     205             : }
     206             : 
     207             : inline
     208         611 : const std::vector<double>& GridVessel::getGridSpacing() const {
     209         611 :   if( gtype==flat ) return dx;
     210           0 :   plumed_merror("dont understand what spacing means for spherical grids");
     211             : }
     212             : 
     213             : inline
     214           4 : double GridVessel::getCellVolume() const {
     215           4 :   if( gtype==flat ) {
     216           8 :     double myvol=1.0; for(unsigned i=0; i<dimension; ++i) myvol *= dx[i];
     217             :     return myvol;
     218             :   } else {
     219           1 :     return 4*pi / static_cast<double>( getNumberOfPoints() );
     220             :   }
     221             : }
     222             : 
     223             : inline
     224             : unsigned GridVessel::getDimension() const {
     225      938454 :   return dimension;
     226             : }
     227             : 
     228             : inline
     229             : bool GridVessel::isPeriodic( const unsigned& i ) const {
     230             :   plumed_dbg_assert( gtype==flat );
     231        7859 :   return pbc[i];
     232             : }
     233             : 
     234             : inline
     235             : std::string GridVessel::getComponentName( const unsigned& i ) const {
     236     1012634 :   return arg_names[i];
     237             : }
     238             : 
     239             : inline
     240          47 : unsigned GridVessel::getNumberOfComponents() const {
     241          47 :   if( noderiv ) return nper;
     242          32 :   return nper / ( dimension + 1 );
     243             : }
     244             : 
     245             : inline
     246             : double GridVessel::getGridExtent( const unsigned& i ) const {
     247             :   plumed_dbg_assert( gtype==flat );
     248          40 :   return max[i] - min[i];
     249             : }
     250             : 
     251             : inline
     252             : bool GridVessel::noDerivatives() const {
     253       14312 :   return noderiv;
     254             : }
     255             : 
     256             : inline
     257    21151408 : bool GridVessel::inactive( const unsigned& ip ) const {
     258             :   plumed_dbg_assert( ip<npoints );
     259    21151408 :   return !active[ip];
     260             : }
     261             : 
     262             : inline
     263             : const std::vector<unsigned>& GridVessel::getStride() const {
     264             :   plumed_dbg_assert( gtype==flat );
     265             :   return stride;
     266             : }
     267             : 
     268             : inline
     269          71 : unsigned GridVessel::getNumberOfBufferPoints() const {
     270         321 :   return npoints;
     271             : }
     272             : 
     273             : inline
     274       26251 : std::string GridVessel::getType() const {
     275       26251 :   if( gtype==flat ) return "flat";
     276         422 :   else if( gtype==fibonacci ) return "fibonacci";
     277           0 :   plumed_error();
     278             : }
     279             : 
     280             : inline
     281           2 : double GridVessel::getFibonacciCutoff() const {
     282           2 :   return 0.0;
     283             : }
     284             : 
     285             : }
     286             : }
     287             : #endif