Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "VesselRegister.h"
23 : #include "Vessel.h"
24 : #include "StoreDataVessel.h"
25 : #include "ActionWithVessel.h"
26 :
27 : namespace PLMD {
28 : namespace vesselbase {
29 :
30 : // This is not the most efficient implementation
31 : // The calculation of all the colvars is parallelized
32 : // but the loops for calculating moments are not
33 : // Feel free to reimplement this if you know how
34 : class Moments : public Vessel {
35 : private:
36 : unsigned mycomponent;
37 : StoreDataVessel* mystash;
38 : std::vector<unsigned> powers;
39 : std::vector<Value*> value_out;
40 : public:
41 : static void registerKeywords( Keywords& keys );
42 : static void reserveKeyword( Keywords& keys );
43 : explicit Moments( const vesselbase::VesselOptions& da );
44 : std::string description() override;
45 : void resize() override;
46 6936 : void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const override {}
47 : void finish( const std::vector<double>& buffer ) override;
48 : bool applyForce( std::vector<double>& forces ) override;
49 : };
50 :
51 13793 : PLUMED_REGISTER_VESSEL(Moments,"MOMENTS")
52 :
53 8 : void Moments::registerKeywords( Keywords& keys ) {
54 8 : Vessel::registerKeywords( keys );
55 8 : keys.remove("LABEL");
56 16 : keys.add("compulsory","COMPONENT","1","the component of the vector for which to calculate this quantity");
57 16 : keys.add("optional","MOMENTS","the list of moments that you would like to calculate");
58 8 : }
59 :
60 4595 : void Moments::reserveKeyword( Keywords& keys ) {
61 9190 : keys.reserve("optional","MOMENTS","calculate the moments of the distribution of collective variables. "
62 : "The mth 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 "
63 : "the average for the distribution. The moments keyword takes a lists of integers as input or a range. Each integer is a value of \\f$m\\f$. The final "
64 : "calculated values can be referenced using moment-\\f$m\\f$. You can use the COMPONENT keyword in this action but the syntax is slightly different. "
65 : "If you would like the second and third moments of the third component you would use MOMENTS={COMPONENT=3 MOMENTS=2-3}. The moments would then be referred to "
66 : "using the labels moment-3-2 and moment-3-3. This syntax is also required if you are using numbered MOMENT keywords i.e. MOMENTS1, MOMENTS2...");
67 9190 : keys.reset_style("MOMENTS","vessel");
68 9190 : keys.addOutputComponent("moment","MOMENTS","the central moments of the distribution of values. The second moment "
69 : "would be referenced elsewhere in the input file using "
70 : "<em>label</em>.moment-2, the third as <em>label</em>.moment-3, etc.");
71 4595 : }
72 :
73 8 : Moments::Moments( const vesselbase::VesselOptions& da) :
74 8 : Vessel(da) {
75 8 : mystash = getAction()->buildDataStashes( NULL );
76 8 : ActionWithValue* a=dynamic_cast<ActionWithValue*>( getAction() );
77 8 : plumed_massert(a,"cannot create passable values as base action does not inherit from ActionWithValue");
78 :
79 : std::vector<std::string> moments;
80 8 : std::string valstr = "moment-";
81 8 : parse("COMPONENT",mycomponent);
82 16 : parseVector("MOMENTS",moments);
83 8 : if( getNumericalLabel()!=0 ) {
84 : std::string numstr;
85 2 : Tools::convert( mycomponent, numstr);
86 4 : valstr = "moment-" + numstr + "-";
87 : }
88 8 : if( moments.size()==0 ) {
89 12 : moments=Tools::getWords(getAllInput(),"\t\n ,");
90 : }
91 8 : Tools::interpretRanges(moments);
92 : unsigned nn;
93 20 : for(unsigned i=0; i<moments.size(); ++i) {
94 24 : a->addComponentWithDerivatives( valstr + moments[i] );
95 24 : a->componentIsNotPeriodic( valstr + moments[i] );
96 12 : value_out.push_back( a->copyOutput( a->getNumberOfComponents()-1 ) );
97 12 : Tools::convert( moments[i], nn );
98 12 : if( nn<2 ) {
99 0 : error("moments are only possible for m>=2" );
100 : }
101 12 : powers.push_back( nn );
102 : std::string num;
103 12 : Tools::convert(powers[i],num);
104 : }
105 8 : }
106 :
107 28 : void Moments::resize() {
108 : resizeBuffer(0);
109 28 : if( getAction()->derivativesAreRequired() ) {
110 56 : for(unsigned i=0; i<value_out.size(); ++i) {
111 36 : value_out[i]->resizeDerivatives( getAction()->getNumberOfDerivatives() );
112 : }
113 : }
114 28 : }
115 :
116 8 : std::string Moments::description() {
117 : std::string descri, num;
118 8 : Tools::convert(powers[0],num);
119 8 : if( getNumericalLabel()==0 ) {
120 12 : descri = "value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
121 8 : for(unsigned i=1; i<powers.size(); ++i) {
122 2 : Tools::convert(powers[i],num);
123 4 : descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
124 : }
125 : } else {
126 : std::string numlab;
127 2 : Tools::convert( mycomponent, numlab );
128 4 : descri = "value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
129 4 : for(unsigned i=1; i<powers.size(); ++i) {
130 2 : Tools::convert(powers[i],num);
131 4 : descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
132 : }
133 : }
134 8 : return descri;
135 : }
136 :
137 822 : void Moments::finish( const std::vector<double>& buffer ) {
138 : const double pi=3.141592653589793238462643383279502884197169399375105820974944592307;
139 : unsigned nvals=getAction()->getFullNumberOfTasks();
140 822 : std::vector<double> myvalues( getAction()->getNumberOfQuantities() );
141 :
142 : double mean=0;
143 822 : Value myvalue;
144 822 : if( getAction()->isPeriodic() ) {
145 : std::string str_min, str_max;
146 0 : getAction()->retrieveDomain( str_min, str_max );
147 : double pfactor, min, max;
148 0 : Tools::convert(str_min,min);
149 0 : Tools::convert(str_max,max);
150 0 : pfactor = 2*pi / ( max-min );
151 0 : myvalue.setDomain( str_min, str_max );
152 : double sinsum=0, cossum=0, val;
153 0 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
154 0 : mystash->retrieveSequentialValue( i, false, myvalues );
155 0 : val=pfactor*( myvalues[mycomponent] - min );
156 0 : sinsum+=std::sin(val);
157 0 : cossum+=std::cos(val);
158 : }
159 0 : mean = 0.5 + std::atan2( sinsum / static_cast<double>( nvals ), cossum / static_cast<double>( nvals ) ) / (2*pi);
160 0 : mean = min + (max-min)*mean;
161 : } else {
162 7758 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
163 6936 : mystash->retrieveSequentialValue( i, false, myvalues );
164 6936 : mean+=myvalues[mycomponent];
165 : }
166 822 : mean/=static_cast<double>( nvals );
167 822 : myvalue.setNotPeriodic();
168 : }
169 :
170 2056 : for(unsigned npow=0; npow<powers.size(); ++npow) {
171 : double dev1=0;
172 1234 : if( value_out[0]->getNumberOfDerivatives()>0 ) {
173 4606 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
174 4172 : mystash->retrieveSequentialValue( i, false, myvalues );
175 4172 : dev1+=pow( myvalue.difference( mean, myvalues[mycomponent] ), powers[npow] - 1 );
176 : }
177 434 : dev1/=static_cast<double>( nvals );
178 : }
179 :
180 : double moment=0;
181 1234 : MultiValue myvals( getAction()->getNumberOfQuantities(), getAction()->getNumberOfDerivatives() );
182 1234 : myvals.clearAll();
183 11006 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
184 9772 : mystash->retrieveSequentialValue( i, false, myvalues );
185 9772 : double tmp=myvalue.difference( mean, myvalues[mycomponent] );
186 9772 : moment+=pow( tmp, powers[npow] );
187 9772 : if( value_out[npow]->getNumberOfDerivatives() ) {
188 4172 : double pref=pow( tmp, powers[npow] - 1 ) - dev1;
189 4172 : mystash->retrieveDerivatives( i, false, myvals );
190 166880 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
191 : unsigned jatom=myvals.getActiveIndex(j);
192 162708 : value_out[npow]->addDerivative(jatom, pref*myvals.getDerivative( mycomponent, jatom ) );
193 : }
194 4172 : myvals.clearAll();
195 : }
196 : }
197 1234 : if( value_out[npow]->getNumberOfDerivatives()>0 ) {
198 434 : value_out[npow]->chainRule( powers[npow] / static_cast<double>( nvals ) );
199 : }
200 1234 : value_out[npow]->set( moment / static_cast<double>( nvals ) );
201 1234 : }
202 1644 : }
203 :
204 432 : bool Moments::applyForce( std::vector<double>& forces ) {
205 432 : std::vector<double> tmpforce( forces.size() );
206 432 : forces.assign(forces.size(),0.0);
207 : bool wasforced=false;
208 1276 : for(unsigned i=0; i<value_out.size(); ++i) {
209 844 : if( value_out[i]->applyForce( tmpforce ) ) {
210 : wasforced=true;
211 0 : for(unsigned j=0; j<forces.size(); ++j) {
212 0 : forces[j]+=tmpforce[j];
213 : }
214 : }
215 : }
216 432 : return wasforced;
217 : }
218 :
219 : }
220 : }
|
