Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2020 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_FunctionOfGrid_h
23 : #define __PLUMED_gridtools_FunctionOfGrid_h
24 :
25 : #include "ActionWithGrid.h"
26 : #include "function/Custom.h"
27 : #include "tools/Matrix.h"
28 :
29 : namespace PLMD {
30 : namespace gridtools {
31 :
32 : template <class T>
33 : class FunctionOfGrid : public ActionWithGrid {
34 : private:
35 : /// The function that is being computed
36 : T myfunc;
37 : public:
38 : static void registerKeywords(Keywords&);
39 : explicit FunctionOfGrid(const ActionOptions&);
40 : /// This does setup required on first step
41 : void setupOnFirstStep( const bool incalc ) override ;
42 : /// Get the number of derivatives for this action
43 : unsigned getNumberOfDerivatives() override ;
44 : /// Get the label to write in the graph
45 0 : std::string writeInGraph() const override {
46 0 : return myfunc.getGraphInfo( getName() );
47 : }
48 : /// Get the underlying names
49 : std::vector<std::string> getGridCoordinateNames() const override ;
50 : /// Get the underlying grid coordinates object
51 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
52 : /// Calculate the function
53 : void performTask( const unsigned& current, MultiValue& myvals ) const override ;
54 : ///
55 : void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
56 : const unsigned& bufstart, std::vector<double>& buffer ) const override ;
57 : /// Add the forces
58 : void apply() override;
59 : };
60 :
61 : template <class T>
62 1047 : void FunctionOfGrid<T>::registerKeywords(Keywords& keys ) {
63 1047 : ActionWithGrid::registerKeywords(keys);
64 1047 : std::string name = keys.getDisplayName();
65 1047 : std::size_t und=name.find("_GRID");
66 2094 : keys.setDisplayName( name.substr(0,und) );
67 1047 : keys.reserve("compulsory","PERIODIC","if the output of your function is periodic then you should specify the periodicity of the function. If the output is not periodic you must state this using PERIODIC=NO");
68 871 : T tfunc;
69 1047 : tfunc.registerKeywords( keys );
70 1047 : if( typeid(tfunc)==typeid(function::Custom()) ) {
71 0 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
72 : }
73 2094 : if( keys.getDisplayName()=="INTEGRATE") {
74 294 : keys.setValueDescription("scalar","the numerical integral of the input function over its whole domain");
75 294 : keys.addInputKeyword("compulsory","ARG","grid","the label of the function on a grid that is being integrated");
76 1800 : } else if( keys.getDisplayName()=="SUM") {
77 58 : keys.setValueDescription("scalar","the sum of the value of the function over all the grid points where it has been evaluated");
78 58 : keys.addInputKeyword("compulsory","ARG","grid","the label of the function on a grid from which we are computing a sum");
79 1742 : } else if( keys.outputComponentExists(".#!value") ) {
80 1742 : keys.setValueDescription("grid","the grid obtained by doing an element-wise application of " + keys.getOutputComponentDescription(".#!value") + " to the input grid");
81 1742 : keys.addInputKeyword("compulsory","ARG","scalar/grid","the labels of the scalars and functions on a grid that we are using to compute the required function");
82 : }
83 1918 : }
84 :
85 : template <class T>
86 523 : FunctionOfGrid<T>::FunctionOfGrid(const ActionOptions&ao):
87 : Action(ao),
88 523 : ActionWithGrid(ao) {
89 523 : if( getNumberOfArguments()==0 ) {
90 0 : error("found no arguments");
91 : }
92 : // This will require a fix
93 523 : if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) {
94 0 : error("first input to this action must be a grid");
95 : }
96 : // Get the shape of the input grid
97 523 : std::vector<unsigned> shape( getPntrToArgument(0)->getShape() );
98 937 : for(unsigned i=1; i<getNumberOfArguments(); ++i ) {
99 414 : if( getPntrToArgument(i)->getRank()==0 ) {
100 122 : continue;
101 : }
102 292 : std::vector<unsigned> s( getPntrToArgument(i)->getShape() );
103 292 : if( s.size()!=shape.size() ) {
104 0 : error("mismatch between dimensionalities of input grids");
105 : }
106 : }
107 : // Read the input and do some checks
108 523 : myfunc.read( this );
109 : // Check we are not calculating an integral
110 523 : if( myfunc.zeroRank() ) {
111 90 : shape.resize(0);
112 : }
113 : // Check that derivatives are available
114 90 : if( !myfunc.derivativesImplemented() ) {
115 : error("derivatives have not been implemended for " + getName() );
116 : }
117 : // Get the names of the components
118 523 : std::vector<std::string> components( keywords.getOutputComponents() );
119 : // Create the values to hold the output
120 1046 : if( components.size()!=1 || components[0]!=".#!value" ) {
121 0 : error("functions of grid should only output one grid");
122 : }
123 523 : addValueWithDerivatives( shape );
124 : // Set the periodicities of the output components
125 523 : myfunc.setPeriodicityForOutputs( this );
126 : // Check if we can turn off the derivatives when they are zero
127 433 : if( myfunc.getDerivativeZeroIfValueIsZero() ) {
128 492 : for(int i=0; i<getNumberOfComponents(); ++i) {
129 246 : getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
130 : }
131 : }
132 523 : setupOnFirstStep( false );
133 1046 : }
134 :
135 : template <class T>
136 1045 : void FunctionOfGrid<T>::setupOnFirstStep( const bool incalc ) {
137 : double volume = 1.0;
138 1045 : const GridCoordinatesObject& mygrid = getGridCoordinatesObject();
139 1045 : unsigned npoints = getPntrToArgument(0)->getNumberOfValues();
140 2090 : if( mygrid.getGridType()=="flat" ) {
141 999 : std::vector<unsigned> shape( getGridCoordinatesObject().getNbin(true) );
142 1794 : for(unsigned i=1; i<getNumberOfArguments(); ++i ) {
143 795 : if( getPntrToArgument(i)->getRank()==0 ) {
144 222 : continue;
145 : }
146 573 : std::vector<unsigned> s( getPntrToArgument(i)->getShape() );
147 1160 : for(unsigned j=0; j<shape.size(); ++j) {
148 587 : if( shape[j]!=s[j] ) {
149 0 : error("mismatch between sizes of input grids");
150 : }
151 : }
152 : }
153 1998 : for(int i=0; i<getNumberOfComponents(); ++i) {
154 999 : if( getPntrToComponent(i)->getRank()>0 ) {
155 833 : getPntrToComponent(i)->setShape(shape);
156 : }
157 : }
158 999 : std::vector<double> vv( getGridCoordinatesObject().getGridSpacing() );
159 166 : volume=vv[0];
160 1081 : for(unsigned i=1; i<vv.size(); ++i) {
161 4 : volume *=vv[i];
162 : }
163 : } else {
164 14 : volume=4*pi / static_cast<double>( npoints );
165 : }
166 : // This resizes the scalars
167 2090 : for(int i=0; i<getNumberOfComponents(); ++i) {
168 1045 : if( getPntrToComponent(i)->getRank()==0 ) {
169 180 : getPntrToComponent(i)->resizeDerivatives( npoints );
170 : }
171 : }
172 1045 : if( getName()=="SUM_GRID" ) {
173 : volume = 1.0;
174 : }
175 : // This sets the prefactor to the volume which converts integrals to sums
176 1045 : myfunc.setup( this );
177 180 : myfunc.setPrefactor( this, volume );
178 1045 : }
179 :
180 : template <class T>
181 15460 : const GridCoordinatesObject& FunctionOfGrid<T>::getGridCoordinatesObject() const {
182 15460 : ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
183 15460 : plumed_assert( ag );
184 15460 : return ag->getGridCoordinatesObject();
185 : }
186 :
187 : template <class T>
188 198 : std::vector<std::string> FunctionOfGrid<T>::getGridCoordinateNames() const {
189 198 : ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
190 198 : plumed_assert( ag );
191 198 : return ag->getGridCoordinateNames();
192 : }
193 :
194 : template <class T>
195 1814 : unsigned FunctionOfGrid<T>::getNumberOfDerivatives() {
196 514 : if( myfunc.zeroRank() ) {
197 514 : return getPntrToArgument(0)->getNumberOfValues();
198 : }
199 1300 : unsigned nder = getGridCoordinatesObject().getDimension();
200 1300 : return getGridCoordinatesObject().getDimension() + getNumberOfArguments() - myfunc.getArgStart();
201 : }
202 :
203 : template <class T>
204 974814 : void FunctionOfGrid<T>::performTask( const unsigned& current, MultiValue& myvals ) const {
205 : unsigned argstart=myfunc.getArgStart();
206 974814 : std::vector<double> args( getNumberOfArguments() - argstart );
207 2629182 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
208 1654368 : if( getPntrToArgument(i)->getRank()==0 ) {
209 544232 : args[i-argstart]=getPntrToArgument(i)->get();
210 : } else {
211 1110136 : args[i-argstart] = getPntrToArgument(i)->get(current);
212 : }
213 : }
214 : // Calculate the function and its derivatives
215 974814 : std::vector<double> vals(1);
216 974814 : Matrix<double> derivatives( 1, getNumberOfArguments()-argstart );
217 974814 : myfunc.calc( this, args, vals, derivatives );
218 974814 : unsigned np = myvals.getTaskIndex();
219 : // And set the values and derivatives
220 974814 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
221 974814 : myvals.addValue( ostrn, vals[0] );
222 23665 : if( !myfunc.zeroRank() ) {
223 : // Add the derivatives for a grid
224 2581852 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
225 : // We store all the derivatives of all the input values - i.e. the grid points these are used in apply
226 1630703 : myvals.addDerivative( ostrn, getConstPntrToComponent(0)->getRank()+j-argstart, derivatives(0,j-argstart) );
227 : // And now we calculate the derivatives of the value that is stored on the grid correctly so that we can interpolate functions
228 1630703 : if( getPntrToArgument(j)->getRank()!=0 ) {
229 3413677 : for(unsigned k=0; k<getPntrToArgument(j)->getRank(); ++k) {
230 2327206 : myvals.addDerivative( ostrn, k, derivatives(0,j-argstart)*getPntrToArgument(j)->getGridDerivative( np, k ) );
231 : }
232 : }
233 : }
234 951149 : unsigned nderivatives = getConstPntrToComponent(0)->getNumberOfGridDerivatives();
235 4575808 : for(unsigned j=0; j<nderivatives; ++j) {
236 3624659 : myvals.updateIndex( ostrn, j );
237 : }
238 23665 : } else if( !doNotCalculateDerivatives() ) {
239 : // These are the derivatives of the integral
240 8161 : myvals.addDerivative( ostrn, current, derivatives(0,0) );
241 8161 : myvals.updateIndex( ostrn, current );
242 : }
243 974814 : }
244 :
245 : template <class T>
246 951149 : void FunctionOfGrid<T>::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
247 : const unsigned& bufstart, std::vector<double>& buffer ) const {
248 951149 : if( getConstPntrToComponent(0)->getRank()>0 && getConstPntrToComponent(0)->hasDerivatives() ) {
249 : plumed_dbg_assert( getNumberOfComponents()==1 && valindex==0 );
250 951149 : unsigned nder = getConstPntrToComponent(0)->getNumberOfGridDerivatives();
251 951149 : unsigned ostr = getConstPntrToComponent(0)->getPositionInStream();
252 951149 : unsigned kp = bufstart + code*(1+nder);
253 951149 : buffer[kp] += myvals.get( ostr );
254 4575808 : for(unsigned i=0; i<nder; ++i) {
255 3624659 : buffer[kp + 1 + i] += myvals.getDerivative( ostr, i );
256 : }
257 : } else {
258 0 : ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer );
259 : }
260 951149 : }
261 :
262 : template <class T>
263 6067 : void FunctionOfGrid<T>::apply() {
264 6067 : if( doNotCalculateDerivatives() || !getPntrToComponent(0)->forcesWereAdded() ) {
265 5442 : return;
266 : }
267 :
268 : // This applies forces for the integral
269 494 : if( myfunc.zeroRank() ) {
270 494 : ActionWithVector::apply();
271 494 : return;
272 : }
273 :
274 : // Work out how to deal with arguments
275 : unsigned nscalars=0, argstart=myfunc.getArgStart();
276 1870 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
277 1245 : if( getPntrToArgument(i)->getRank()==0 ) {
278 20 : nscalars++;
279 : }
280 : }
281 :
282 625 : std::vector<double> totv(nscalars,0);
283 625 : Value* outval=getPntrToComponent(0);
284 7575 : for(unsigned i=0; i<outval->getNumberOfValues(); ++i) {
285 : nscalars=0;
286 19845 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
287 12895 : double fforce = outval->getForce(i);
288 12895 : if( getPntrToArgument(j)->getRank()==0 ) {
289 3205 : totv[nscalars] += fforce*outval->getGridDerivative( i, outval->getRank()+j );
290 3205 : nscalars++;
291 : } else {
292 9690 : double vval = outval->getGridDerivative( i, outval->getRank()+j );
293 9690 : getPntrToArgument(j)->addForce( i, fforce*vval );
294 : }
295 : }
296 : }
297 : nscalars=0;
298 1870 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
299 1245 : if( getPntrToArgument(i)->getRank()==0 ) {
300 20 : getPntrToArgument(i)->addForce( 0, totv[nscalars] );
301 20 : nscalars++;
302 : }
303 : }
304 :
305 : }
306 :
307 : }
308 : }
309 : #endif
|
