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 : {
79 : bool normalize;
80 : std::vector<double> coefficients;
81 : std::vector<double> parameters;
82 : std::vector<double> powers;
83 : public:
84 : explicit Combine(const ActionOptions&);
85 : void calculate() override;
86 : static void registerKeywords(Keywords& keys);
87 : };
88 :
89 :
90 10690 : PLUMED_REGISTER_ACTION(Combine,"COMBINE")
91 :
92 138 : void Combine::registerKeywords(Keywords& keys) {
93 138 : Function::registerKeywords(keys);
94 276 : keys.use("ARG"); keys.use("PERIODIC");
95 276 : keys.add("compulsory","COEFFICIENTS","1.0","the coefficients of the arguments in your function");
96 276 : keys.add("compulsory","PARAMETERS","0.0","the parameters of the arguments in your function");
97 276 : keys.add("compulsory","POWERS","1.0","the powers to which you are raising each of the arguments in your function");
98 276 : keys.addFlag("NORMALIZE",false,"normalize all the coefficients so that in total they are equal to one");
99 138 : }
100 :
101 137 : Combine::Combine(const ActionOptions&ao):
102 : Action(ao),
103 : Function(ao),
104 137 : normalize(false),
105 277 : coefficients(getNumberOfArguments(),1.0),
106 140 : parameters(getNumberOfArguments(),0.0),
107 277 : powers(getNumberOfArguments(),1.0)
108 : {
109 273 : parseVector("COEFFICIENTS",coefficients);
110 136 : if(coefficients.size()!=static_cast<unsigned>(getNumberOfArguments()))
111 0 : error("Size of COEFFICIENTS array should be the same as number for arguments");
112 :
113 271 : parseVector("PARAMETERS",parameters);
114 135 : if(parameters.size()!=static_cast<unsigned>(getNumberOfArguments()))
115 0 : error("Size of PARAMETERS array should be the same as number for arguments");
116 :
117 269 : parseVector("POWERS",powers);
118 134 : if(powers.size()!=static_cast<unsigned>(getNumberOfArguments()))
119 0 : error("Size of POWERS array should be the same as number for arguments");
120 :
121 137 : parseFlag("NORMALIZE",normalize);
122 :
123 134 : if(normalize) {
124 : double n=0.0;
125 16 : for(unsigned i=0; i<coefficients.size(); i++) n+=coefficients[i];
126 16 : for(unsigned i=0; i<coefficients.size(); i++) coefficients[i]*=(1.0/n);
127 : }
128 :
129 134 : addValueWithDerivatives();
130 134 : checkRead();
131 :
132 134 : log.printf(" with coefficients:");
133 404 : for(unsigned i=0; i<coefficients.size(); i++) log.printf(" %f",coefficients[i]);
134 134 : log.printf("\n");
135 134 : log.printf(" with parameters:");
136 404 : for(unsigned i=0; i<parameters.size(); i++) log.printf(" %f",parameters[i]);
137 134 : log.printf("\n");
138 134 : log.printf(" and powers:");
139 404 : for(unsigned i=0; i<powers.size(); i++) log.printf(" %f",powers[i]);
140 134 : log.printf("\n");
141 140 : }
142 :
143 7887 : void Combine::calculate() {
144 7887 : double combine=0.0;
145 20242 : for(unsigned i=0; i<coefficients.size(); ++i) {
146 12355 : double cv = (getArgument(i)-parameters[i]);
147 12355 : combine+=coefficients[i]*std::pow(cv,powers[i]);
148 12355 : setDerivative(i,coefficients[i]*powers[i]*std::pow(cv,powers[i]-1.0));
149 : };
150 7887 : setValue(combine);
151 7887 : }
152 :
153 : }
154 : }
155 :
156 :
|
