Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 : #ifndef __PLUMED_function_FunctionSetup_h
23 : #define __PLUMED_function_FunctionSetup_h
24 :
25 : #include "core/ActionWithValue.h"
26 : #include "core/ActionWithArguments.h"
27 : #include "core/ActionWithVector.h"
28 :
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace function {
33 :
34 4527 : struct FunctionOptions {
35 : /// Is the derivative zero if the value is zero
36 : bool derivativeZeroIfValueIsZero = false;
37 : /// Are multiple components registered with a single name as in SphericalHarmonic and Moments
38 : std::vector<std::string> multipleValuesForEachRegisteredComponent;
39 : };
40 :
41 : template <class T>
42 2232 : struct FunctionData {
43 : /// Set equal to one if we are doing EvaluateGridFunction
44 : unsigned argstart = 0;
45 : // Number of scalars that appear in input
46 : unsigned nscalars = 0;
47 : T f;
48 : // This is for setting up the functions
49 : static void setup( T& myfunc,
50 : const std::vector<std::string>& components,
51 : const std::vector<std::size_t>& shape,
52 : bool hasderiv,
53 : ActionWithValue* action );
54 : #ifdef __PLUMED_HAS_OPENACC
55 : void toACCDevice() const {
56 : #pragma acc enter data copyin(this[0:1], argstart, nscalars)
57 : f.toACCDevice();
58 : }
59 : void removeFromACCDevice() const {
60 : f.removeFromACCDevice();
61 : #pragma acc exit data delete(nscalars, argstart, this[0:1])
62 : }
63 : #endif //__PLUMED_HAS_OPENACC
64 : };
65 :
66 : template <class T>
67 4395 : void FunctionData<T>::setup( T& myfunc,
68 : const std::vector<std::string>& components,
69 : const std::vector<std::size_t>& shape,
70 : bool hasderiv,
71 : ActionWithValue* action ) {
72 4395 : ActionWithArguments* aarg = action->castToActionWithArguments();
73 4395 : plumed_assert( aarg );
74 4395 : FunctionOptions options;
75 4395 : T::read( myfunc, aarg, options );
76 :
77 4388 : if( action->keywords.getDisplayName()=="SORT" ) {
78 41 : for(unsigned j=0; j<aarg->getNumberOfArguments(); ++j) {
79 : std::string num;
80 28 : Tools::convert( j+1, num );
81 28 : if( hasderiv ) {
82 26 : action->addComponentWithDerivatives( num, shape );
83 : } else {
84 2 : action->addComponent( num, shape );
85 : }
86 : }
87 : } else {
88 8796 : for(unsigned i=0; i<components.size(); ++i) {
89 4421 : if( options.multipleValuesForEachRegisteredComponent.size()>0 ) {
90 840 : for(unsigned j=0; j<options.multipleValuesForEachRegisteredComponent.size(); ++j) {
91 743 : if( hasderiv ) {
92 30 : action->addComponentWithDerivatives( components[i] + options.multipleValuesForEachRegisteredComponent[j], shape );
93 : } else {
94 1456 : action->addComponent( components[i] + options.multipleValuesForEachRegisteredComponent[j], shape );
95 : }
96 : }
97 4324 : } else if( components[i]==".#!value" ) {
98 4318 : if( hasderiv ) {
99 2300 : action->addValueWithDerivatives( shape );
100 : } else {
101 2018 : action->addValue( shape );
102 : }
103 6 : } else if( components[i].find_first_of("_")!=std::string::npos ) {
104 2 : if( aarg->getNumberOfArguments()==1 && hasderiv ) {
105 1 : action->addValueWithDerivatives( shape );
106 1 : } else if( aarg->getNumberOfArguments()==1 ) {
107 0 : action->addValue( shape );
108 : } else {
109 3 : for(unsigned j=0; j<aarg->getNumberOfArguments(); ++j) {
110 2 : if( hasderiv ) {
111 4 : action->addComponentWithDerivatives( (aarg->getPntrToArgument(j))->getName() + components[i], shape );
112 : } else {
113 0 : action->addComponent( (aarg->getPntrToArgument(j))->getName() + components[i], shape );
114 : }
115 : }
116 : }
117 : } else {
118 4 : if( hasderiv ) {
119 0 : action->addComponentWithDerivatives( components[i], shape );
120 : } else {
121 4 : action->addComponent( components[i], shape );
122 : }
123 : }
124 : }
125 : }
126 4388 : if( options.derivativeZeroIfValueIsZero ) {
127 3143 : for(unsigned i=0; i<action->getNumberOfComponents(); ++i) {
128 1914 : (action->copyOutput(i))->setDerivativeIsZeroWhenValueIsZero();
129 : }
130 : }
131 4388 : if( action->keywords.exists("PERIODIC") ) {
132 3766 : plumed_assert( action->keywords.getDisplayName()!="DIFFERENCE" );
133 3766 : std::vector<std::string> period;
134 7532 : action->parseVector("PERIODIC",period);
135 3766 : if( period.size()==1 ) {
136 3758 : if( period[0]!="NO") {
137 0 : action->error("input to PERIODIC keyword does not make sense");
138 : }
139 7516 : for(unsigned i=0; i<action->getNumberOfComponents(); ++i) {
140 3758 : (action->copyOutput(i))->setNotPeriodic();
141 : }
142 : return;
143 8 : } else if( period.size()!=2 ) {
144 0 : action->error("input to PERIODIC keyword does not make sense");
145 : }
146 16 : for(unsigned i=0; i<action->getNumberOfComponents(); ++i) {
147 8 : (action->copyOutput(i))->setDomain( period[0], period[1] );
148 : }
149 4395 : } else if( action->keywords.getDisplayName()=="DIFFERENCE" ) {
150 184 : if( aarg->getPntrToArgument(0)->isPeriodic() ) {
151 : std::string min, max;
152 42 : aarg->getPntrToArgument(0)->getDomain( min, max );
153 42 : action->setPeriodic(min,max);
154 142 : } else if( aarg->getPntrToArgument(1)->isPeriodic() ) {
155 : std::string min, max;
156 0 : aarg->getPntrToArgument(1)->getDomain( min, max );
157 0 : action->setPeriodic( min, max );
158 : } else {
159 142 : action->setNotPeriodic();
160 : }
161 : } else {
162 1584 : for(unsigned i=0; i<action->getNumberOfComponents(); ++i) {
163 1146 : (action->copyOutput(i))->setNotPeriodic();
164 : }
165 : }
166 : }
167 :
168 : struct FunctionOutput {
169 : View<double> values;
170 : View2D<double> derivs;
171 : static FunctionOutput create( unsigned nvals,
172 : double* vals,
173 : unsigned ndev,
174 : double* devs ) {
175 : return FunctionOutput{
176 : View<double>(vals,nvals),
177 : View2D<double>(devs,nvals,ndev)
178 124668 : };
179 : }
180 : };
181 :
182 : /// Check if all the values returned by the function are zero when all the arguments are zero.
183 : /// When this is not the case an element of a vector cannot be skipped because its arguments are zero.
184 : /// A value that is not finite (e.g. 0/0) is treated as zero, as that element is not meaningful.
185 : template <class T>
186 1290 : bool functionIsZeroWhenArgumentsAreZero( const T& f, std::size_t nargs, std::size_t ncomponents ) {
187 1290 : std::vector<double> zeros( nargs, 0.0 ), vals( ncomponents, 0.0 ), derivs( ncomponents*nargs, 0.0 );
188 : auto funcout = FunctionOutput::create( ncomponents, vals.data(), nargs, derivs.data() );
189 1290 : T::calc( f, true, View<const double>( zeros.data(), nargs ), funcout );
190 2369 : for(const auto & v : vals) {
191 1291 : if( std::isfinite(v) && std::fabs(v)>epsilon ) {
192 : return false;
193 : }
194 : }
195 : return true;
196 : }
197 :
198 : } // namespace function
199 : } // namespace PLMD
200 : #endif
|
