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 : #include "Function.h"
23 : #include "ActionRegister.h"
24 :
25 : namespace PLMD {
26 : namespace function {
27 :
28 : //+PLUMEDOC FUNCTION COMBINE
29 : /*
30 : Calculate a polynomial combination of a set of other variables.
31 :
32 : The functional form of this function is
33 : \f[
34 : C=\sum_{i=1}^{N_{arg}} c_i (x_i-a_i)^{p_i}
35 : \f]
36 :
37 : The coefficients c, the parameters a and the powers p are provided as vectors.
38 :
39 : Notice that COMBINE is not able to predict which will be periodic domain
40 : of the computed value automatically. The user is thus forced to specify it
41 : explicitly. Use PERIODIC=NO if the resulting variable is not periodic,
42 : and PERIODIC=A,B where A and B are the two boundaries if the resulting variable
43 : is periodic.
44 :
45 :
46 :
47 : \par Examples
48 :
49 : The following input tells plumed to print the distance between atoms 3 and 5
50 : its square (as computed from the x,y,z components) and the distance
51 : again as computed from the square root of the square.
52 : \plumedfile
53 : DISTANCE LABEL=dist ATOMS=3,5 COMPONENTS
54 : COMBINE LABEL=distance2 ARG=dist.x,dist.y,dist.z POWERS=2,2,2 PERIODIC=NO
55 : COMBINE LABEL=distance ARG=distance2 POWERS=0.5 PERIODIC=NO
56 : PRINT ARG=distance,distance2
57 : \endplumedfile
58 : (See also \ref PRINT and \ref DISTANCE).
59 :
60 : The following input tells plumed to add a restraint on the
61 : cube of a dihedral angle. Notice that since the angle has a
62 : periodic domain
63 : -pi,pi its cube has a domain -pi**3,pi**3.
64 : \plumedfile
65 : t: TORSION ATOMS=1,3,5,7
66 : c: COMBINE ARG=t POWERS=3 PERIODIC=-31.0062766802998,31.0062766802998
67 : RESTRAINT ARG=c KAPPA=10 AT=0
68 : \endplumedfile
69 :
70 :
71 :
72 : */
73 : //+ENDPLUMEDOC
74 :
75 :
76 : class Combine :
77 : public Function {
78 : bool normalize;
79 : std::vector<double> coefficients;
80 : std::vector<double> parameters;
81 : std::vector<double> powers;
82 : public:
83 : explicit Combine(const ActionOptions&);
84 : void calculate() override;
85 : static void registerKeywords(Keywords& keys);
86 : };
87 :
88 :
89 14062 : PLUMED_REGISTER_ACTION(Combine,"COMBINE")
90 :
91 144 : void Combine::registerKeywords(Keywords& keys) {
92 144 : Function::registerKeywords(keys);
93 144 : keys.use("ARG");
94 144 : keys.use("PERIODIC");
95 288 : keys.add("compulsory","COEFFICIENTS","1.0","the coefficients of the arguments in your function");
96 288 : keys.add("compulsory","PARAMETERS","0.0","the parameters of the arguments in your function");
97 288 : keys.add("compulsory","POWERS","1.0","the powers to which you are raising each of the arguments in your function");
98 288 : keys.addFlag("NORMALIZE",false,"normalize all the coefficients so that in total they are equal to one");
99 144 : }
100 :
101 140 : Combine::Combine(const ActionOptions&ao):
102 : Action(ao),
103 : Function(ao),
104 140 : normalize(false),
105 283 : coefficients(getNumberOfArguments(),1.0),
106 143 : parameters(getNumberOfArguments(),0.0),
107 283 : powers(getNumberOfArguments(),1.0) {
108 279 : parseVector("COEFFICIENTS",coefficients);
109 139 : if(coefficients.size()!=static_cast<unsigned>(getNumberOfArguments())) {
110 0 : error("Size of COEFFICIENTS array should be the same as number for arguments");
111 : }
112 :
113 277 : parseVector("PARAMETERS",parameters);
114 138 : if(parameters.size()!=static_cast<unsigned>(getNumberOfArguments())) {
115 0 : error("Size of PARAMETERS array should be the same as number for arguments");
116 : }
117 :
118 275 : parseVector("POWERS",powers);
119 137 : if(powers.size()!=static_cast<unsigned>(getNumberOfArguments())) {
120 0 : error("Size of POWERS array should be the same as number for arguments");
121 : }
122 :
123 140 : parseFlag("NORMALIZE",normalize);
124 :
125 137 : if(normalize) {
126 : double n=0.0;
127 16 : for(unsigned i=0; i<coefficients.size(); i++) {
128 14 : n+=coefficients[i];
129 : }
130 16 : for(unsigned i=0; i<coefficients.size(); i++) {
131 14 : coefficients[i]*=(1.0/n);
132 : }
133 : }
134 :
135 137 : addValueWithDerivatives();
136 137 : checkRead();
137 :
138 137 : log.printf(" with coefficients:");
139 413 : for(unsigned i=0; i<coefficients.size(); i++) {
140 276 : log.printf(" %f",coefficients[i]);
141 : }
142 137 : log.printf("\n");
143 137 : log.printf(" with parameters:");
144 413 : for(unsigned i=0; i<parameters.size(); i++) {
145 276 : log.printf(" %f",parameters[i]);
146 : }
147 137 : log.printf("\n");
148 137 : log.printf(" and powers:");
149 413 : for(unsigned i=0; i<powers.size(); i++) {
150 276 : log.printf(" %f",powers[i]);
151 : }
152 137 : log.printf("\n");
153 143 : }
154 :
155 7905 : void Combine::calculate() {
156 7905 : double combine=0.0;
157 20296 : for(unsigned i=0; i<coefficients.size(); ++i) {
158 12391 : double cv = (getArgument(i)-parameters[i]);
159 12391 : combine+=coefficients[i]*std::pow(cv,powers[i]);
160 12391 : setDerivative(i,coefficients[i]*powers[i]*std::pow(cv,powers[i]-1.0));
161 : };
162 7905 : setValue(combine);
163 7905 : }
164 :
165 : }
166 : }
167 :
168 :
|
