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 161 : 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 161 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
81 161 : keys.addFlag("MIDPOINTS",false,"interpolate the values of the function at the midpoints of the grid coordinates of the input grid");
82 : EvaluateGridFunction ii;
83 161 : ii.registerKeywords( keys );
84 322 : keys.setValueDescription("grid","the function evaluated onto the interpolated grid");
85 161 : }
86 :
87 82 : InterpolateGrid::InterpolateGrid(const ActionOptions&ao):
88 : Action(ao),
89 82 : ActionWithGrid(ao) {
90 82 : if( getNumberOfArguments()!=1 ) {
91 0 : error("should only be one argument to this action");
92 : }
93 82 : if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) {
94 0 : error("input to this action should be a grid");
95 : }
96 :
97 82 : parseFlag("MIDPOINTS",midpoints);
98 82 : parseVector("GRID_BIN",nbin);
99 164 : parseVector("GRID_SPACING",gspacing);
100 : unsigned dimension = getPntrToArgument(0)->getRank();
101 82 : if( !midpoints && nbin.size()!=dimension && gspacing.size()!=dimension ) {
102 0 : error("MIDPOINTS, GRID_BIN or GRID_SPACING must be set");
103 : }
104 82 : if( midpoints ) {
105 80 : log.printf(" evaluating function at midpoints of cells in input grid\n");
106 2 : } else if( nbin.size()==dimension ) {
107 2 : log.printf(" number of bins in grid %d", nbin[0]);
108 2 : for(unsigned i=1; i<nbin.size(); ++i) {
109 0 : log.printf(", %d", nbin[i]);
110 : }
111 2 : log.printf("\n");
112 0 : } else if( gspacing.size()==dimension ) {
113 0 : log.printf(" spacing for bins in grid %f", gspacing[0]);
114 0 : for(unsigned i=1; i<gspacing.size(); ++i) {
115 0 : log.printf(", %d", gspacing[i]);
116 : }
117 0 : log.printf("\n");
118 : }
119 : // Create the input grid
120 82 : input_grid.read( this );
121 : // Need this for creation of tasks
122 164 : output_grid.setup( "flat", input_grid.getPbc(), 0, 0.0 );
123 :
124 : // Now add a value
125 82 : std::vector<unsigned> shape( dimension, 0 );
126 82 : addValueWithDerivatives( shape );
127 :
128 82 : if( getPntrToArgument(0)->isPeriodic() ) {
129 : std::string min, max;
130 0 : getPntrToArgument(0)->getDomain( min, max );
131 0 : setPeriodic( min, max );
132 : } else {
133 82 : setNotPeriodic();
134 : }
135 82 : setupOnFirstStep( false );
136 82 : }
137 :
138 164 : void InterpolateGrid::setupOnFirstStep( const bool incalc ) {
139 164 : input_grid.setup( this );
140 164 : ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
141 164 : plumed_assert( ag );
142 164 : const GridCoordinatesObject& mygrid = ag->getGridCoordinatesObject();
143 164 : if( midpoints ) {
144 : double min, max;
145 160 : nbin.resize( getPntrToComponent(0)->getRank() );
146 : std::vector<std::string> str_min( input_grid.getMin() ), str_max(input_grid.getMax() );
147 320 : for(unsigned i=0; i<nbin.size(); ++i) {
148 160 : if( incalc ) {
149 80 : Tools::convert( str_min[i], min );
150 80 : Tools::convert( str_max[i], max );
151 80 : min += 0.5*input_grid.getGridSpacing()[i];
152 : }
153 320 : if( input_grid.getPbc()[i] ) {
154 30 : nbin[i] = input_grid.getNbin()[i];
155 30 : if( incalc ) {
156 15 : max += 0.5*input_grid.getGridSpacing()[i];
157 : }
158 : } else {
159 130 : nbin[i] = input_grid.getNbin()[i] - 1;
160 130 : if( incalc ) {
161 65 : max -= 0.5*input_grid.getGridSpacing()[i];
162 : }
163 : }
164 160 : if( incalc ) {
165 80 : Tools::convert( min, str_min[i] );
166 80 : Tools::convert( max, str_max[i] );
167 : }
168 : }
169 160 : output_grid.setBounds( str_min, str_max, nbin, gspacing );
170 160 : } else {
171 4 : output_grid.setBounds( mygrid.getMin(), mygrid.getMax(), nbin, gspacing );
172 : }
173 164 : getPntrToComponent(0)->setShape( output_grid.getNbin(true) );
174 164 : if( !incalc ) {
175 82 : gspacing.resize(0);
176 : }
177 164 : }
178 :
179 256 : unsigned InterpolateGrid::getNumberOfDerivatives() {
180 256 : return getPntrToArgument(0)->getRank();
181 : }
182 :
183 392 : const GridCoordinatesObject& InterpolateGrid::getGridCoordinatesObject() const {
184 392 : return output_grid;
185 : }
186 :
187 2 : std::vector<std::string> InterpolateGrid::getGridCoordinateNames() const {
188 2 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
189 2 : plumed_assert( ag );
190 2 : return ag->getGridCoordinateNames();
191 : }
192 :
193 15616 : void InterpolateGrid::performTask( const unsigned& current, MultiValue& myvals ) const {
194 15616 : std::vector<double> pos( output_grid.getDimension() );
195 15616 : output_grid.getGridPointCoordinates( current, pos );
196 15616 : std::vector<double> val(1);
197 : Matrix<double> der( 1, output_grid.getDimension() );
198 15616 : input_grid.calc( this, pos, val, der );
199 15616 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
200 15616 : myvals.setValue( ostrn, val[0] );
201 31232 : for(unsigned i=0; i<output_grid.getDimension(); ++i) {
202 15616 : myvals.addDerivative( ostrn, i, der(0,i) );
203 15616 : myvals.updateIndex( ostrn, i );
204 : }
205 15616 : }
206 :
207 13660 : void InterpolateGrid::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
208 : const unsigned& bufstart, std::vector<double>& buffer ) const {
209 : plumed_dbg_assert( valindex==0 );
210 13660 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
211 13660 : unsigned istart = bufstart + (1+output_grid.getDimension())*code;
212 13660 : buffer[istart] += myvals.get( ostrn );
213 27320 : for(unsigned i=0; i<output_grid.getDimension(); ++i) {
214 13660 : buffer[istart+1+i] += myvals.getDerivative( ostrn, i );
215 : }
216 13660 : }
217 :
218 1956 : void InterpolateGrid::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
219 1956 : std::vector<double> pos(output_grid.getDimension());
220 1956 : double ff = myval->getForce(itask);
221 1956 : output_grid.getGridPointCoordinates( itask, pos );
222 1956 : input_grid.applyForce( this, pos, ff, forces );
223 1956 : }
224 :
225 :
226 : }
227 : }
|
