Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 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 24 : 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();
45 : void resize();
46 6936 : void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {}
47 : void finish( const std::vector<double>& buffer );
48 : bool applyForce( std::vector<double>& forces );
49 : };
50 :
51 7372 : PLUMED_REGISTER_VESSEL(Moments,"MOMENTS")
52 :
53 8 : void Moments::registerKeywords( Keywords& keys ) {
54 8 : Vessel::registerKeywords( keys );
55 16 : keys.remove("LABEL");
56 40 : keys.add("compulsory","COMPONENT","1","the component of the vector for which to calculate this quantity");
57 32 : keys.add("compulsory","MOMENTS","the list of moments that you would like to calculate");
58 8 : }
59 :
60 1839 : void Moments::reserveKeyword( Keywords& keys ) {
61 7356 : 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 5517 : keys.reset_style("MOMENTS","vessel");
68 7356 : 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 1839 : }
72 :
73 8 : Moments::Moments( const vesselbase::VesselOptions& da) :
74 16 : Vessel(da)
75 : {
76 8 : mystash = getAction()->buildDataStashes( NULL );
77 8 : ActionWithValue* a=dynamic_cast<ActionWithValue*>( getAction() );
78 8 : plumed_massert(a,"cannot create passable values as base action does not inherit from ActionWithValue");
79 :
80 8 : std::vector<std::string> moments; std::string valstr;
81 8 : if( getNumericalLabel()==0 ) {
82 : valstr = "moment-";
83 18 : moments=Tools::getWords(getAllInput(),"\t\n ,");
84 6 : Tools::interpretRanges(moments); mycomponent=1;
85 : } else {
86 4 : std::string numstr; parse("COMPONENT",mycomponent);
87 6 : Tools::convert( mycomponent, numstr); valstr = "moment-" + numstr + "-";
88 4 : parseVector("MOMENTS",moments); Tools::interpretRanges(moments);
89 : }
90 : unsigned nn;
91 52 : for(unsigned i=0; i<moments.size(); ++i) {
92 24 : a->addComponentWithDerivatives( valstr + moments[i] );
93 24 : a->componentIsNotPeriodic( valstr + moments[i] );
94 24 : value_out.push_back( a->copyOutput( a->getNumberOfComponents()-1 ) );
95 12 : Tools::convert( moments[i], nn );
96 12 : if( nn<2 ) error("moments are only possible for m>=2" );
97 24 : powers.push_back( nn ); std::string num; Tools::convert(powers[i],num);
98 : }
99 8 : }
100 :
101 28 : void Moments::resize() {
102 : resizeBuffer(0);
103 28 : if( getAction()->derivativesAreRequired() ) {
104 220 : for(unsigned i=0; i<value_out.size(); ++i) value_out[i]->resizeDerivatives( getAction()->getNumberOfDerivatives() );
105 : }
106 28 : }
107 :
108 8 : std::string Moments::description() {
109 : std::string descri, num;
110 8 : Tools::convert(powers[0],num);
111 8 : if( getNumericalLabel()==0 ) {
112 54 : descri = "value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
113 18 : for(unsigned i=1; i<powers.size(); ++i) {
114 2 : Tools::convert(powers[i],num);
115 20 : descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
116 : }
117 : } else {
118 2 : std::string numlab; Tools::convert( mycomponent, numlab );
119 24 : descri = "value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
120 10 : for(unsigned i=1; i<powers.size(); ++i) {
121 2 : Tools::convert(powers[i],num);
122 26 : descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
123 : }
124 : }
125 8 : return descri;
126 : }
127 :
128 822 : void Moments::finish( const std::vector<double>& buffer ) {
129 : const double pi=3.141592653589793238462643383279502884197169399375105820974944592307;
130 : unsigned nvals=getAction()->getFullNumberOfTasks();
131 822 : std::vector<double> myvalues( getAction()->getNumberOfQuantities() );
132 :
133 1644 : double mean=0; Value myvalue;
134 822 : if( getAction()->isPeriodic() ) {
135 0 : std::string str_min, str_max; getAction()->retrieveDomain( str_min, str_max );
136 0 : double pfactor, min, max; Tools::convert(str_min,min); Tools::convert(str_max,max);
137 0 : pfactor = 2*pi / ( max-min ); myvalue.setDomain( str_min, str_max );
138 : double sinsum=0, cossum=0, val;
139 0 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
140 0 : mystash->retrieveSequentialValue( i, false, myvalues );
141 0 : val=pfactor*( myvalues[mycomponent] - min );
142 0 : sinsum+=sin(val); cossum+=cos(val);
143 : }
144 0 : mean = 0.5 + atan2( sinsum / static_cast<double>( nvals ), cossum / static_cast<double>( nvals ) ) / (2*pi);
145 0 : mean = min + (max-min)*mean;
146 : } else {
147 7758 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
148 13872 : mystash->retrieveSequentialValue( i, false, myvalues ); mean+=myvalues[mycomponent];
149 : }
150 822 : mean/=static_cast<double>( nvals ); myvalue.setNotPeriodic();
151 : }
152 :
153 5346 : for(unsigned npow=0; npow<powers.size(); ++npow) {
154 : double dev1=0;
155 1234 : if( value_out[0]->getNumberOfDerivatives()>0 ) {
156 4606 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
157 4172 : mystash->retrieveSequentialValue( i, false, myvalues );
158 8344 : dev1+=pow( myvalue.difference( mean, myvalues[mycomponent] ), powers[npow] - 1 );
159 : }
160 434 : dev1/=static_cast<double>( nvals );
161 : }
162 :
163 : double moment=0;
164 3702 : MultiValue myvals( getAction()->getNumberOfQuantities(), getAction()->getNumberOfDerivatives() ); myvals.clearAll();
165 11006 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
166 9772 : mystash->retrieveSequentialValue( i, false, myvalues );
167 19544 : double tmp=myvalue.difference( mean, myvalues[mycomponent] );
168 9772 : moment+=pow( tmp, powers[npow] );
169 9772 : if( value_out[npow]->getNumberOfDerivatives() ) {
170 4172 : double pref=pow( tmp, powers[npow] - 1 ) - dev1;
171 4172 : mystash->retrieveDerivatives( i, false, myvals );
172 329588 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
173 : unsigned jatom=myvals.getActiveIndex(j);
174 325416 : value_out[npow]->addDerivative(jatom, pref*myvals.getDerivative( mycomponent, jatom ) );
175 : }
176 4172 : myvals.clearAll();
177 : }
178 : }
179 2102 : if( value_out[npow]->getNumberOfDerivatives()>0 ) value_out[npow]->chainRule( powers[npow] / static_cast<double>( nvals ) );
180 1234 : value_out[npow]->set( moment / static_cast<double>( nvals ) );
181 : }
182 822 : }
183 :
184 432 : bool Moments::applyForce( std::vector<double>& forces ) {
185 432 : std::vector<double> tmpforce( forces.size() );
186 864 : forces.assign(forces.size(),0.0); bool wasforced=false;
187 3396 : for(unsigned i=0; i<value_out.size(); ++i) {
188 844 : if( value_out[i]->applyForce( tmpforce ) ) {
189 : wasforced=true;
190 0 : for(unsigned j=0; j<forces.size(); ++j) forces[j]+=tmpforce[j];
191 : }
192 : }
193 432 : return wasforced;
194 : }
195 :
196 : }
197 5517 : }
|
