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 "FunctionWithSingleArgument.h"
26 : #include "core/ActionRegister.h"
27 : #include "FunctionSetup.h"
28 :
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace function {
33 :
34 : //+PLUMEDOC FUNCTION MOMENTS
35 : /*
36 : Calculate central moments from the distribution of input quantities
37 :
38 : 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:
39 :
40 : $$
41 : \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
42 : $$
43 :
44 : 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
45 : and third central moment for the distribution of the four input distances.
46 :
47 : ```plumed
48 : d12: DISTANCE ATOMS=1,2
49 : d13: DISTANCE ATOMS=1,3
50 : d14: DISTANCE ATOMS=1,4
51 : d15: DISTANCE ATOMS=1,5
52 : mv: MOMENTS ARG=d12,d13,d14,d15 POWERS=2,3
53 : PRINT ARG=mv.moment-2,mv.moment-3 FILE=colvar
54 : ```
55 :
56 : Notice that you can also achieve the same result using the following input:
57 :
58 : ```plumed
59 : d: DISTANCE ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5
60 : mv: MOMENTS ARG=d POWERS=2,3
61 : PRINT ARG=mv.moment-2,mv.moment-3 FILE=colvar
62 : ```
63 :
64 : In this second case the four distances are passed to the MOMENTS action as a vector. The MOMENTS action then outputs 2 components - the
65 : two central moments that were requested.
66 :
67 : 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.
68 : It does not accept matrices or a list of vectors in input.
69 :
70 : */
71 : //+ENDPLUMEDOC
72 :
73 1 : class Moments {
74 : public:
75 : bool isperiodic;
76 : double min, max, pfactor;
77 : double max_minus_min;
78 : double inv_max_minus_min;
79 : std::vector<int> powers;
80 : static void registerKeywords(Keywords& keys);
81 : static void read( Moments& func, ActionWithArguments* action, FunctionOptions& options );
82 : static void calc( const Moments& func, bool noderiv, View<const double> args, FunctionOutput& funcout );
83 : };
84 :
85 : typedef FunctionShortcut<Moments> MomentsShortcut;
86 : PLUMED_REGISTER_ACTION(MomentsShortcut,"MOMENTS")
87 : typedef FunctionWithSingleArgument<Moments> SingleargMoments;
88 : PLUMED_REGISTER_ACTION(SingleargMoments,"MOMENTS_ONEARG")
89 : typedef FunctionOfScalar<Moments> ScalarMoments;
90 : PLUMED_REGISTER_ACTION(ScalarMoments,"MOMENTS_SCALAR")
91 : typedef FunctionOfVector<Moments> VectorMoments;
92 : PLUMED_REGISTER_ACTION(VectorMoments,"MOMENTS_VECTOR")
93 :
94 53 : void Moments::registerKeywords(Keywords& keys) {
95 53 : keys.add("compulsory","POWERS","calculate the central moments of the distribution of collective variables. "
96 : "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 "
97 : "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 "
98 : "calculated values can be referenced using moment-\\f$m\\f$.");
99 106 : keys.addOutputComponent("moment","default","scalar","the central moments of the distribution of values. The second central moment "
100 : "would be referenced elsewhere in the input file using "
101 : "<em>label</em>.moment-2, the third as <em>label</em>.moment-3, etc.");
102 53 : }
103 :
104 9 : void Moments::read( Moments& func, ActionWithArguments* action, FunctionOptions& options ) {
105 9 : func.isperiodic = (action->getPntrToArgument(0))->isPeriodic();
106 9 : if( func.isperiodic ) {
107 : std::string str_min, str_max;
108 2 : (action->getPntrToArgument(0))->getDomain( str_min, str_max );
109 2 : for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
110 0 : if( !(action->getPntrToArgument(i))->isPeriodic() ) {
111 0 : action->error("cannot mix periodic and non periodic variables when calculating moments");
112 : }
113 : std::string str_min2, str_max2;
114 0 : (action->getPntrToArgument(i))->getDomain( str_min2, str_max2);
115 0 : if( str_min!=str_min2 || str_max!=str_max2 ) {
116 0 : action->error("all input arguments should have same domain when calculating moments");
117 : }
118 : }
119 2 : Tools::convert(str_min,func.min);
120 2 : Tools::convert(str_max,func.max);
121 2 : func.max_minus_min = func.max - func.min;
122 2 : func.inv_max_minus_min = 1 / func.max_minus_min;
123 2 : func.pfactor = 2*pi*func.inv_max_minus_min;
124 : } else {
125 16 : for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
126 9 : if( (action->getPntrToArgument(i))->isPeriodic() ) {
127 0 : action->error("cannot mix periodic and non periodic variables when calculating moments");
128 : }
129 : }
130 : }
131 :
132 9 : action->parseVector("POWERS",func.powers);
133 24 : for(unsigned i=0; i<func.powers.size(); ++i) {
134 15 : if( func.powers[i]<2 ) {
135 0 : action->error("first central moment is zero do you really need to calculate that");
136 : }
137 15 : action->log.printf(" computing %dth central moment of distribution of input cvs \n", func.powers[i]);
138 : std::string num;
139 15 : Tools::convert(func.powers[i],num);
140 30 : options.multipleValuesForEachRegisteredComponent.push_back( "-" + num );
141 : }
142 9 : }
143 :
144 234 : void Moments::calc( const Moments& func, bool noderiv, const View<const double> args, FunctionOutput& funcout ) {
145 : double mean=0;
146 234 : double inorm = 1.0 / static_cast<double>( args.size() );
147 234 : if( func.isperiodic ) {
148 : // Circular mean: the angles are measured from the lower end of the domain,
149 : // so the mean is min + (max-min)*atan2(<sin>,<cos>)/(2 pi)
150 : double sinsum=0, cossum=0;
151 60 : for(unsigned i=0; i<args.size(); ++i) {
152 50 : const double val=func.pfactor*( args[i] - func.min );
153 50 : sinsum+=sin(val);
154 50 : cossum+=cos(val);
155 : }
156 10 : const double meanangle = atan2( sinsum, cossum );
157 10 : mean = func.min + func.max_minus_min*meanangle/(2*pi);
158 : // derivative of the circular mean with respect to argument i is cos(theta_i-theta_mean)/R
159 : // with R=|sum_j exp(i theta_j)|; it reduces to 1/N when all the angles are equal
160 10 : const double r2 = sinsum*sinsum + cossum*cossum;
161 10 : const double invr = r2>0 ? 1.0/std::sqrt(r2) : 0.0;
162 30 : for(unsigned npow=0; npow<func.powers.size(); ++npow) {
163 : double dev1=0;
164 120 : for(unsigned i=0; i<args.size(); ++i) {
165 100 : double s = func.inv_max_minus_min*( args[i] - mean );
166 100 : s = Tools::pbc(s);
167 100 : dev1+=pow( s*func.max_minus_min, func.powers[npow] - 1 );
168 : }
169 20 : funcout.values[npow] = 0;
170 20 : double prefactor = func.powers[npow]*inorm;
171 120 : for(unsigned i=0; i<args.size(); ++i) {
172 100 : double tmp = func.inv_max_minus_min*( args[i] - mean );
173 100 : tmp = func.max_minus_min*Tools::pbc(tmp);
174 100 : funcout.values[npow] += inorm*pow( tmp, func.powers[npow] );
175 100 : if( !noderiv ) {
176 100 : const double dmean = invr*cos( func.pfactor*( args[i] - func.min ) - meanangle );
177 100 : funcout.derivs[npow][i] = prefactor*(pow( tmp, func.powers[npow] - 1 ) - dev1*dmean);
178 : }
179 : }
180 : }
181 : } else {
182 1768 : for(unsigned i=0; i<args.size(); ++i) {
183 1544 : mean+=args[i];
184 : }
185 224 : mean *= inorm;
186 661 : for(unsigned npow=0; npow<func.powers.size(); ++npow) {
187 : double dev1=0;
188 3421 : for(unsigned i=0; i<args.size(); ++i) {
189 2984 : dev1+=pow( args[i] - mean, func.powers[npow] - 1 );
190 : }
191 437 : dev1*=inorm;
192 437 : funcout.values[npow] = 0;
193 437 : double prefactor = func.powers[npow]*inorm;
194 3421 : for(unsigned i=0; i<args.size(); ++i) {
195 2984 : double tmp=args[i] - mean;
196 2984 : funcout.values[npow] += inorm*pow( tmp, func.powers[npow] );
197 2984 : if( !noderiv ) {
198 112 : funcout.derivs[npow][i] = prefactor*(pow( tmp, func.powers[npow] - 1 ) - dev1);
199 : }
200 : }
201 : }
202 : }
203 234 : }
204 :
205 : }
206 : }
207 :
208 :
|
