Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "core/ActionRegister.h"
23 : #include "core/PlumedMain.h"
24 : #include "EvaluateGridFunction.h"
25 : #include "ActionWithGrid.h"
26 :
27 : //+PLUMEDOC GRIDANALYSIS INTERPOLATE_GRID
28 : /*
29 : Interpolate a smooth function stored on a grid onto a grid with a smaller grid spacing.
30 :
31 : This action takes a function evaluated on a grid as input and can be used to interpolate the values of that
32 : function on to a finer grained grid. The interpolation within this algorithm is done using splines.
33 :
34 : \par Examples
35 :
36 : The input below can be used to post process a trajectory. It calculates a \ref HISTOGRAM as a function the
37 : distance between atoms 1 and 2 using kernel density estimation. During the calculation the values of the kernels
38 : are evaluated at 100 points on a uniform grid between 0.0 and 3.0. Prior to outputting this function at the end of the
39 : simulation this function is interpolated onto a finer grid of 200 points between 0.0 and 3.0.
40 :
41 : \plumedfile
42 : x: DISTANCE ATOMS=1,2
43 : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
44 : ii: INTERPOLATE_GRID GRID=hA1 GRID_BIN=200
45 : DUMPGRID GRID=ii FILE=histo.dat
46 : \endplumedfile
47 :
48 : */
49 : //+ENDPLUMEDOC
50 :
51 : namespace PLMD {
52 : namespace gridtools {
53 :
54 : class InterpolateGrid : public ActionWithGrid {
55 : private:
56 : bool midpoints;
57 : std::vector<unsigned> nbin;
58 : std::vector<double> gspacing;
59 : EvaluateGridFunction input_grid;
60 : GridCoordinatesObject output_grid;
61 : public:
62 : static void registerKeywords( Keywords& keys );
63 : explicit InterpolateGrid(const ActionOptions&ao);
64 : void setupOnFirstStep( const bool incalc ) override ;
65 : unsigned getNumberOfDerivatives() override ;
66 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
67 : std::vector<std::string> getGridCoordinateNames() const override ;
68 : void performTask( const unsigned& current, MultiValue& myvals ) const override ;
69 : void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
70 : const unsigned& bufstart, std::vector<double>& buffer ) const ;
71 : void gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const override ;
72 : };
73 :
74 : PLUMED_REGISTER_ACTION(InterpolateGrid,"INTERPOLATE_GRID")
75 :
76 161 : void InterpolateGrid::registerKeywords( Keywords& keys ) {
77 161 : ActionWithGrid::registerKeywords( keys );
78 322 : keys.add("optional","GRID_BIN","the number of bins for the grid");
79 322 : keys.addInputKeyword("compulsory","ARG","grid","the label for function on the grid that you would like to interpolate");
80 322 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
81 322 : keys.addFlag("MIDPOINTS",false,"interpolate the values of the function at the midpoints of the grid coordinates of the input grid");
82 161 : EvaluateGridFunction ii; ii.registerKeywords( keys );
83 322 : keys.setValueDescription("grid","the function evaluated onto the interpolated grid");
84 161 : }
85 :
86 82 : InterpolateGrid::InterpolateGrid(const ActionOptions&ao):
87 : Action(ao),
88 82 : ActionWithGrid(ao)
89 : {
90 82 : if( getNumberOfArguments()!=1 ) error("should only be one argument to this action");
91 82 : if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) error("input to this action should be a grid");
92 :
93 328 : parseFlag("MIDPOINTS",midpoints); parseVector("GRID_BIN",nbin); parseVector("GRID_SPACING",gspacing); unsigned dimension = getPntrToArgument(0)->getRank();
94 82 : if( !midpoints && nbin.size()!=dimension && gspacing.size()!=dimension ) error("MIDPOINTS, GRID_BIN or GRID_SPACING must be set");
95 82 : if( midpoints ) {
96 80 : log.printf(" evaluating function at midpoints of cells in input grid\n");
97 2 : } else if( nbin.size()==dimension ) {
98 2 : log.printf(" number of bins in grid %d", nbin[0]);
99 2 : for(unsigned i=1; i<nbin.size(); ++i) log.printf(", %d", nbin[i]);
100 2 : log.printf("\n");
101 0 : } else if( gspacing.size()==dimension ) {
102 0 : log.printf(" spacing for bins in grid %f", gspacing[0]);
103 0 : for(unsigned i=1; i<gspacing.size(); ++i) log.printf(", %d", gspacing[i]);
104 0 : log.printf("\n");
105 : }
106 : // Create the input grid
107 82 : input_grid.read( this );
108 : // Need this for creation of tasks
109 164 : output_grid.setup( "flat", input_grid.getPbc(), 0, 0.0 );
110 :
111 : // Now add a value
112 82 : std::vector<unsigned> shape( dimension, 0 );
113 82 : addValueWithDerivatives( shape );
114 :
115 82 : if( getPntrToArgument(0)->isPeriodic() ) {
116 0 : std::string min, max; getPntrToArgument(0)->getDomain( min, max ); setPeriodic( min, max );
117 82 : } else setNotPeriodic();
118 82 : setupOnFirstStep( false );
119 82 : }
120 :
121 164 : void InterpolateGrid::setupOnFirstStep( const bool incalc ) {
122 164 : input_grid.setup( this );
123 164 : ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
124 164 : plumed_assert( ag ); const GridCoordinatesObject& mygrid = ag->getGridCoordinatesObject();
125 164 : if( midpoints ) {
126 160 : double min, max; nbin.resize( getPntrToComponent(0)->getRank() );
127 : std::vector<std::string> str_min( input_grid.getMin() ), str_max(input_grid.getMax() );
128 320 : for(unsigned i=0; i<nbin.size(); ++i) {
129 160 : if( incalc ) {
130 80 : Tools::convert( str_min[i], min ); Tools::convert( str_max[i], max );
131 80 : min += 0.5*input_grid.getGridSpacing()[i];
132 : }
133 320 : if( input_grid.getPbc()[i] ) {
134 30 : nbin[i] = input_grid.getNbin()[i];
135 45 : if( incalc ) max += 0.5*input_grid.getGridSpacing()[i];
136 : } else {
137 130 : nbin[i] = input_grid.getNbin()[i] - 1;
138 195 : if( incalc ) max -= 0.5*input_grid.getGridSpacing()[i];
139 : }
140 160 : if( incalc ) {
141 80 : Tools::convert( min, str_min[i] ); Tools::convert( max, str_max[i] );
142 : }
143 : }
144 160 : output_grid.setBounds( str_min, str_max, nbin, gspacing );
145 164 : } else output_grid.setBounds( mygrid.getMin(), mygrid.getMax(), nbin, gspacing );
146 164 : getPntrToComponent(0)->setShape( output_grid.getNbin(true) );
147 164 : if( !incalc ) gspacing.resize(0);
148 164 : }
149 :
150 256 : unsigned InterpolateGrid::getNumberOfDerivatives() {
151 256 : return getPntrToArgument(0)->getRank();
152 : }
153 :
154 392 : const GridCoordinatesObject& InterpolateGrid::getGridCoordinatesObject() const {
155 392 : return output_grid;
156 : }
157 :
158 2 : std::vector<std::string> InterpolateGrid::getGridCoordinateNames() const {
159 2 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
160 2 : plumed_assert( ag ); return ag->getGridCoordinateNames();
161 : }
162 :
163 15616 : void InterpolateGrid::performTask( const unsigned& current, MultiValue& myvals ) const {
164 15616 : std::vector<double> pos( output_grid.getDimension() ); output_grid.getGridPointCoordinates( current, pos );
165 15616 : std::vector<double> val(1); Matrix<double> der( 1, output_grid.getDimension() ); input_grid.calc( this, pos, val, der );
166 15616 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(); myvals.setValue( ostrn, val[0] );
167 31232 : for(unsigned i=0; i<output_grid.getDimension(); ++i) { myvals.addDerivative( ostrn, i, der(0,i) ); myvals.updateIndex( ostrn, i ); }
168 15616 : }
169 :
170 13660 : void InterpolateGrid::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
171 : const unsigned& bufstart, std::vector<double>& buffer ) const {
172 13660 : plumed_dbg_assert( valindex==0 ); unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
173 13660 : unsigned istart = bufstart + (1+output_grid.getDimension())*code; buffer[istart] += myvals.get( ostrn );
174 27320 : for(unsigned i=0; i<output_grid.getDimension(); ++i) buffer[istart+1+i] += myvals.getDerivative( ostrn, i );
175 13660 : }
176 :
177 1956 : void InterpolateGrid::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
178 1956 : std::vector<double> pos(output_grid.getDimension()); double ff = myval->getForce(itask);
179 1956 : output_grid.getGridPointCoordinates( itask, pos ); input_grid.applyForce( this, pos, ff, forces );
180 1956 : }
181 :
182 :
183 : }
184 : }
|
