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