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