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_Combine_h
23 : #define __PLUMED_function_Combine_h
24 : #include "FunctionSetup.h"
25 : #include "tools/Keywords.h"
26 : #include "tools/View.h"
27 : #include "core/ActionWithArguments.h"
28 :
29 : #include <cmath>
30 : #include <numeric>
31 : namespace PLMD {
32 : namespace function {
33 :
34 : class Combine {
35 : struct component {
36 : double coefficient{0.0};
37 : double parameter{0.0};
38 : double power{0.0};
39 : double max_minus_min{0.0};
40 : double inv_max_minus_min{0.0};
41 : bool periodic{false};
42 : /// power is exactly one: avoid the calls to pow
43 : bool linear{false};
44 : };
45 : std::vector<component> components{};
46 : public:
47 : std::size_t ncomponents{0};
48 : component * cmps=nullptr;
49 : /// with CIRCULAR the arguments are combined as angles: atan2(sum c_i sin(p_i theta_i),sum c_i cos(p_i theta_i))
50 : bool circular{false};
51 : /// (max-min)/(2 pi) for the output domain, used with CIRCULAR
52 : double output_scale{1.0};
53 : void update() {
54 492 : cmps = components.data();
55 492 : ncomponents = components.size();
56 : }
57 990 : Combine() = default;
58 369 : ~Combine() = default;
59 : Combine(const Combine&x):
60 : components(x.components),
61 : circular(x.circular),
62 : output_scale(x.output_scale) {
63 : update();
64 : }
65 : Combine(Combine&&x):
66 : components(std::move(x.components)),
67 : circular(x.circular),
68 : output_scale(x.output_scale) {
69 : update();
70 : }
71 : Combine &operator=(const Combine&x) {
72 : if (this!=&x) {
73 : components=x.components;
74 : circular=x.circular;
75 : output_scale=x.output_scale;
76 : update();
77 : }
78 : return *this;
79 : }
80 : Combine &operator=(Combine&&x) {
81 : if (this!=&x) {
82 : components=std::move(x.components);
83 : circular=x.circular;
84 : output_scale=x.output_scale;
85 : update();
86 : }
87 : return *this;
88 : }
89 : static void registerKeywords(Keywords& keys);
90 : static void read( Combine& func,
91 : ActionWithArguments* action,
92 : FunctionOptions& options );
93 : static void calc( const Combine& func,
94 : bool noderiv,
95 : const View<const double>& args,
96 : FunctionOutput& funcout );
97 : #ifdef __PLUMED_HAS_OPENACC
98 : void toACCDevice() const {
99 : #pragma acc enter data copyin(this[0:1], ncomponents, cmps[0:ncomponents])
100 : }
101 : void removeFromACCDevice() const {
102 : #pragma acc exit data delete(cmps[0:ncomponents], ncomponents, this[0:1])
103 : }
104 : #endif // __PLUMED_HAS_OPENACC
105 : };
106 :
107 1808 : void Combine::registerKeywords(Keywords& keys) {
108 1808 : keys.use("PERIODIC");
109 1808 : keys.add("compulsory","COEFFICIENTS","1.0","the coefficients of the arguments in your function");
110 1808 : keys.add("compulsory","PARAMETERS","0.0","the parameters of the arguments in your function");
111 1808 : keys.add("compulsory","POWERS","1.0","the powers to which you are raising each of the arguments in your function");
112 1808 : keys.addFlag("NORMALIZE",false,"normalize all the coefficients so that in total they are equal to one");
113 1808 : keys.addFlag("CIRCULAR",false,"all the arguments are periodic and are combined as angles: the output is the argument of the complex number "
114 : "sum_i c_i exp(i p_i theta_i), where theta_i is (x_i-a_i) converted to an angle, scaled to the domain given with PERIODIC. "
115 : "With the default coefficients and powers this is the circular mean of the arguments");
116 3616 : keys.setValueDescription("scalar/vector/matrix","a linear combination");
117 1808 : }
118 :
119 495 : void Combine::read( Combine& func, ActionWithArguments* action,
120 : FunctionOptions& options ) {
121 495 : unsigned nargs = action->getNumberOfArguments();
122 495 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(action);
123 495 : if(av && av->getNumberOfMasks()>0) {
124 18 : nargs = nargs - av->getNumberOfMasks();
125 : }
126 495 : std::vector<double> coefficients(nargs);
127 989 : action->parseVector("COEFFICIENTS",coefficients);
128 :
129 494 : if(coefficients.size()!=nargs) {
130 0 : action->error("Size of COEFFICIENTS array should be the same as number for arguments");
131 : }
132 497 : std::vector<double> parameters(nargs);
133 987 : action->parseVector("PARAMETERS",parameters);
134 493 : if(parameters.size()!=nargs) {
135 0 : action->error("Size of PARAMETERS array should be the same as number for arguments");
136 : }
137 495 : std::vector<double> powers(nargs);
138 985 : action->parseVector("POWERS",powers);
139 492 : if(powers.size()!=nargs) {
140 0 : action->error("Size of POWERS array should be the same as number for arguments");
141 : }
142 :
143 : bool normalize;
144 492 : action->parseFlag("NORMALIZE",normalize);
145 492 : if(normalize) {
146 : double n=std::accumulate(coefficients.begin(),coefficients.end(),0.0);
147 : std::transform(coefficients.begin(),coefficients.end(),coefficients.begin(),
148 : [n](double x) {
149 38 : return x/n;
150 : });
151 : }
152 :
153 492 : action->log.printf(" with coefficients:");
154 492 : func.components.resize( nargs );
155 : func.update();
156 2005 : for(unsigned i=0; i<func.components.size(); i++) {
157 1513 : func.components[i].coefficient = coefficients[i];
158 1513 : func.components[i].parameter = parameters[i];
159 1513 : func.components[i].power = powers[i];
160 1513 : action->log.printf(" %f",func.components[i].coefficient);
161 : }
162 492 : action->log.printf("\n with parameters:");
163 2005 : for(unsigned i=0; i<func.components.size(); i++) {
164 1513 : action->log.printf(" %f",func.components[i].parameter);
165 : }
166 492 : action->log.printf("\n and powers:");
167 2005 : for(unsigned i=0; i<func.components.size(); i++) {
168 1513 : action->log.printf(" %f",func.components[i].power);
169 : }
170 492 : action->log.printf("\n");
171 2005 : for(unsigned i=0; i<nargs; ++i) {
172 1513 : func.components[i].linear = (func.components[i].power==1.0);
173 : }
174 : // Store periodicity stuff
175 2005 : for(unsigned i=0; i<nargs; ++i) {
176 1513 : if( (action->getPntrToArgument(i))->isPeriodic() ) {
177 62 : func.components[i].periodic = true;
178 : std::string min, max;
179 62 : (action->getPntrToArgument(i))->getDomain( min, max );
180 : double dmin, dmax;
181 62 : Tools::convert( min, dmin );
182 62 : Tools::convert( max, dmax );
183 62 : func.components[i].max_minus_min = dmax - dmin;
184 62 : func.components[i].inv_max_minus_min = 1.0 / func.components[i].max_minus_min ;
185 : }
186 : }
187 :
188 492 : action->parseFlag("CIRCULAR",func.circular);
189 492 : if( func.circular ) {
190 11 : for(unsigned i=0; i<nargs; ++i) {
191 8 : if( !func.components[i].periodic ) {
192 0 : action->error("with CIRCULAR all the arguments must be periodic (" + action->getPntrToArgument(i)->getName() + " is not)");
193 : }
194 8 : const double p=func.components[i].power;
195 8 : if( p<0 || std::fabs(p-std::round(p))>0.0 ) {
196 0 : action->error("with CIRCULAR the POWERS must be non-negative integers");
197 : }
198 : }
199 : // The output domain is needed to scale the angle; PERIODIC is read later by FunctionData::setup
200 4 : std::string pstr=action->getKeyword("PERIODIC");
201 6 : std::vector<std::string> period=Tools::getWords(pstr.substr(pstr.find("=")+1),"\t\n ,");
202 3 : if( period.size()!=2 ) {
203 0 : action->error("with CIRCULAR the output is periodic: use PERIODIC=min,max to give its domain");
204 : }
205 : double omin, omax;
206 3 : Tools::convert( period[0], omin );
207 3 : Tools::convert( period[1], omax );
208 3 : if( !(omax>omin) ) {
209 0 : action->error("with CIRCULAR the domain given with PERIODIC must have max>min");
210 : }
211 3 : func.output_scale = (omax-omin)/(2*pi);
212 3 : action->log.printf(" combining the arguments as angles (CIRCULAR); output domain %s\n", pstr.c_str());
213 3 : }
214 492 : }
215 :
216 3582321 : void Combine::calc( const Combine& func,
217 : bool noderiv,
218 : const View<const double>& args,
219 : FunctionOutput& funcout ) {
220 3582321 : if( func.circular ) {
221 : // theta_i = 2 pi pbc(x_i-a_i)/(max_i-min_i), output = scale*atan2( sum c_i sin(p_i theta_i), sum c_i cos(p_i theta_i) )
222 : double ssum=0.0, csum=0.0;
223 133 : for(unsigned i=0; i<func.ncomponents; ++i) {
224 92 : const auto& cmp = func.cmps[i];
225 92 : const double angle = cmp.power*2*pi*Tools::pbc( (args[i] - cmp.parameter)*cmp.inv_max_minus_min );
226 92 : ssum += cmp.coefficient*sin(angle);
227 92 : csum += cmp.coefficient*cos(angle);
228 : }
229 41 : funcout.values[0] = func.output_scale*atan2( ssum, csum );
230 41 : if( !noderiv ) {
231 25 : const double r2 = ssum*ssum + csum*csum;
232 : // the result is undefined when the sum vanishes: use zero derivatives there
233 25 : const double ir2 = r2>0 ? 1.0/r2 : 0.0;
234 85 : for(unsigned i=0; i<func.ncomponents; ++i) {
235 60 : const auto& cmp = func.cmps[i];
236 60 : const double angle = cmp.power*2*pi*Tools::pbc( (args[i] - cmp.parameter)*cmp.inv_max_minus_min );
237 60 : funcout.derivs[0][i] = func.output_scale*cmp.coefficient*cmp.power*2*pi*cmp.inv_max_minus_min*
238 60 : ( csum*cos(angle) + ssum*sin(angle) )*ir2;
239 : }
240 : }
241 41 : return;
242 : }
243 3582280 : funcout.values[0]=0.0;
244 14378714 : for(unsigned i=0; i<func.ncomponents; ++i) {
245 10796434 : const auto& cmp = func.cmps[i];
246 10796434 : double cv = args[i] - cmp.parameter;
247 10796434 : if( cmp.periodic ) {
248 2078249 : cv = cmp.max_minus_min*Tools::pbc( cv*cmp.inv_max_minus_min );
249 : }
250 10796434 : if( cmp.linear ) {
251 : // same result as pow(cv,1) and pow(cv,0), without the calls
252 2276964 : funcout.values[0] += cmp.coefficient*cv;
253 2276964 : if( !noderiv ) {
254 73599 : funcout.derivs[0][i] = cmp.coefficient;
255 : }
256 2276964 : continue;
257 : }
258 8519470 : funcout.values[0] += cmp.coefficient*pow( cv, cmp.power );
259 8519470 : if( !noderiv ) {
260 1258209 : funcout.derivs[0][i] = cmp.coefficient*cmp.power*pow(cv,cmp.power-1.0);
261 : }
262 : }
263 : }
264 :
265 : } // namespace function
266 : } // namespace PLMD
267 : #endif //__PLUMED_function_Combine_h
|
