Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2017 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 "FunctionShortcut.h"
23 : #include "FunctionOfScalar.h"
24 : #include "FunctionOfVector.h"
25 : #include "core/ActionRegister.h"
26 : #include "FunctionTemplateBase.h"
27 :
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 : namespace function {
32 :
33 : //+PLUMEDOC FUNCTION MOMENTS
34 : /*
35 : Calculate central moments from the distribution of input quantities
36 :
37 : This action takes a set of $N$ input arguments, $s_i$, and evaluates the $k$th central moment of the distribution of input arguments using:
38 :
39 : $$
40 : \mu_k = \frac{1}{N} \sum_{i=1}^N ( s_i - \langle s \rangle )^k \qquad \textrm{where} \qquad \langle s \rangle = \frac{1}{N} \sum_{i=1}^N s_i
41 : $$
42 :
43 : A single moments action can evaluate more than one central moment at once so, for example, the input below can be used to calculate the second
44 : and third central moment for the distribution of the four input distances.
45 :
46 : ```plumed
47 : d12: DISTANCE ATOMS=1,2
48 : d13: DISTANCE ATOMS=1,3
49 : d14: DISTANCE ATOMS=1,4
50 : d15: DISTANCE ATOMS=1,5
51 : mv: MOMENTS ARG=d12,d13,d14,d15 POWERS=2,3
52 : PRINT ARG=mv.moment-2,mv.moment-3 FILE=colvar
53 : ```
54 :
55 : Notice that you can also achieve the same result using the following input:
56 :
57 : ```plumed
58 : d: DISTANCE ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5
59 : mv: MOMENTS ARG=d POWERS=2,3
60 : PRINT ARG=mv.moment-2,mv.moment-3 FILE=colvar
61 : ```
62 :
63 : In this second case the four distances are passed to the MOMENTS action as a vector. The MOMENTS action then outputs 2 components - the
64 : two central moments that were requested.
65 :
66 : These examples are representative the only two ways you can use this action. In input it can accept either a list of scalars or a single vector.
67 : It does not accept matrices or a list of vectors in input.
68 :
69 : */
70 : //+ENDPLUMEDOC
71 :
72 : //+PLUMEDOC FUNCTION MOMENTS_SCALAR
73 : /*
74 : Calculate the moments of the distribution of input quantities
75 :
76 : \par Examples
77 :
78 : */
79 : //+ENDPLUMEDOC
80 :
81 : //+PLUMEDOC FUNCTION MOMENTS_VECTOR
82 : /*
83 : Calculate the moments of the distribution of input vectors
84 :
85 : \par Examples
86 :
87 : */
88 : //+ENDPLUMEDOC
89 :
90 69 : class Moments : public FunctionTemplateBase {
91 : bool isperiodic, scalar_out;
92 : double min, max, pfactor;
93 : std::vector<int> powers;
94 : public:
95 : void registerKeywords(Keywords& keys) override;
96 : void read( ActionWithArguments* action ) override;
97 0 : bool zeroRank() const override {
98 7 : return scalar_out;
99 : }
100 0 : bool doWithTasks() const override {
101 221 : return !scalar_out;
102 : }
103 : std::vector<std::string> getComponentsPerLabel() const override ;
104 : void setPeriodicityForOutputs( ActionWithValue* action ) override;
105 : void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
106 : };
107 :
108 : typedef FunctionShortcut<Moments> MomentsShortcut;
109 : PLUMED_REGISTER_ACTION(MomentsShortcut,"MOMENTS")
110 : typedef FunctionOfScalar<Moments> ScalarMoments;
111 : PLUMED_REGISTER_ACTION(ScalarMoments,"MOMENTS_SCALAR")
112 : typedef FunctionOfVector<Moments> VectorMoments;
113 : PLUMED_REGISTER_ACTION(VectorMoments,"MOMENTS_VECTOR")
114 :
115 27 : void Moments::registerKeywords(Keywords& keys) {
116 27 : keys.add("compulsory","POWERS","calculate the central moments of the distribution of collective variables. "
117 : "The \\f$m\\f$th central moment of a distribution is calculated using \\f$\\frac{1}{N} \\sum_{i=1}^N ( s_i - \\overline{s} )^m \\f$, where \\f$\\overline{s}\\f$ is "
118 : "the average for the distribution. The POWERS keyword takes a lists of integers as input or a range. Each integer is a value of \\f$m\\f$. The final "
119 : "calculated values can be referenced using moment-\\f$m\\f$.");
120 54 : keys.addOutputComponent("moment","default","scalar","the central moments of the distribution of values. The second central moment "
121 : "would be referenced elsewhere in the input file using "
122 : "<em>label</em>.moment-2, the third as <em>label</em>.moment-3, etc.");
123 27 : }
124 :
125 5 : void Moments::read( ActionWithArguments* action ) {
126 5 : scalar_out = action->getNumberOfArguments()==1;
127 5 : if( scalar_out && action->getPntrToArgument(0)->getRank()==0 ) {
128 0 : action->error("cannot calculate moments if only given one variable");
129 : }
130 :
131 5 : isperiodic = (action->getPntrToArgument(0))->isPeriodic();
132 5 : if( isperiodic ) {
133 : std::string str_min, str_max;
134 0 : (action->getPntrToArgument(0))->getDomain( str_min, str_max );
135 0 : for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
136 0 : if( !(action->getPntrToArgument(i))->isPeriodic() ) {
137 0 : action->error("cannot mix periodic and non periodic variables when calculating moments");
138 : }
139 : std::string str_min2, str_max2;
140 0 : (action->getPntrToArgument(i))->getDomain( str_min2, str_max2);
141 0 : if( str_min!=str_min2 || str_max!=str_max2 ) {
142 0 : action->error("all input arguments should have same domain when calculating moments");
143 : }
144 : }
145 0 : Tools::convert(str_min,min);
146 0 : Tools::convert(str_max,max);
147 0 : pfactor = 2*pi / ( max-min );
148 : } else {
149 14 : for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
150 9 : if( (action->getPntrToArgument(i))->isPeriodic() ) {
151 0 : action->error("cannot mix periodic and non periodic variables when calculating moments");
152 : }
153 : }
154 : }
155 :
156 5 : parseVector(action,"POWERS",powers);
157 13 : for(unsigned i=0; i<powers.size(); ++i) {
158 8 : if( powers[i]<2 ) {
159 0 : action->error("first central moment is zero do you really need to calculate that");
160 : }
161 8 : action->log.printf(" computing %dth central moment of distribution of input cvs \n", powers[i]);
162 : }
163 5 : }
164 :
165 5 : std::vector<std::string> Moments::getComponentsPerLabel() const {
166 : std::vector<std::string> comp;
167 : std::string num;
168 13 : for(unsigned i=0; i<powers.size(); ++i) {
169 8 : Tools::convert(powers[i],num);
170 16 : comp.push_back( "-" + num );
171 : }
172 5 : return comp;
173 0 : }
174 :
175 5 : void Moments::setPeriodicityForOutputs( ActionWithValue* action ) {
176 13 : for(unsigned i=0; i<powers.size(); ++i) {
177 : std::string num;
178 8 : Tools::convert(powers[i],num);
179 16 : action->componentIsNotPeriodic("moment-" + num);
180 : }
181 5 : }
182 :
183 222 : void Moments::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
184 : double mean=0;
185 222 : double inorm = 1.0 / static_cast<double>( args.size() );
186 222 : if( isperiodic ) {
187 : double sinsum=0, cossum=0, val;
188 0 : for(unsigned i=0; i<args.size(); ++i) {
189 0 : val=pfactor*( args[i] - min );
190 0 : sinsum+=sin(val);
191 0 : cossum+=cos(val);
192 : }
193 0 : mean = 0.5 + atan2( inorm*sinsum, inorm*cossum ) / (2*pi);
194 0 : mean = min + (max-min)*mean;
195 : } else {
196 1758 : for(unsigned i=0; i<args.size(); ++i) {
197 1536 : mean+=args[i];
198 : }
199 222 : mean *= inorm;
200 : }
201 :
202 : Value* arg0 = action->getPntrToArgument(0);
203 656 : for(unsigned npow=0; npow<powers.size(); ++npow) {
204 : double dev1=0;
205 3406 : for(unsigned i=0; i<args.size(); ++i) {
206 2972 : dev1+=pow( arg0->difference( mean, args[i] ), powers[npow] - 1 );
207 : }
208 434 : dev1*=inorm;
209 434 : vals[npow] = 0;
210 434 : double prefactor = powers[npow]*inorm;
211 3406 : for(unsigned i=0; i<args.size(); ++i) {
212 2972 : double tmp=arg0->difference( mean, args[i] );
213 2972 : vals[npow] += inorm*pow( tmp, powers[npow] );
214 2972 : derivatives(npow,i) = prefactor*(pow( tmp, powers[npow] - 1 ) - dev1);
215 : }
216 : }
217 222 : }
218 :
219 : }
220 : }
221 :
222 :
|
