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