Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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 "StoreDataVessel.h"
24 : #include "ActionWithVessel.h"
25 :
26 : namespace PLMD {
27 : namespace vesselbase {
28 :
29 : // This is not the most efficient implementation
30 : // The calculation of all the colvars is parallelized
31 : // but the loops for calculating moments are not
32 : // Feel free to reimplement this if you know how
33 10 : class Moments : public StoreDataVessel {
34 : private:
35 : std::vector<unsigned> powers;
36 : std::vector<Value*> value_out;
37 : public:
38 : static void registerKeywords( Keywords& keys );
39 : static void reserveKeyword( Keywords& keys );
40 : explicit Moments( const vesselbase::VesselOptions& da );
41 : std::string description();
42 : void resize();
43 : void finish( const std::vector<double>& buffer );
44 : bool applyForce( std::vector<double>& forces );
45 : };
46 :
47 2528 : PLUMED_REGISTER_VESSEL(Moments,"MOMENTS")
48 :
49 5 : void Moments::registerKeywords( Keywords& keys ) {
50 5 : StoreDataVessel::registerKeywords( keys );
51 5 : }
52 :
53 841 : void Moments::reserveKeyword( Keywords& keys ) {
54 : keys.reserve("optional","MOMENTS","calculate the moments of the distribution of collective variables. "
55 : "The \\f$m\\f$th 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 "
56 : "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 "
57 841 : "calculated values can be referenced using moment-\\f$m\\f$.");
58 841 : keys.reset_style("MOMENTS","vessel");
59 : keys.addOutputComponent("moment","MOMENTS","the central moments of the distribution of values. The second moment "
60 : "would be referenced elsewhere in the input file using "
61 841 : "<em>label</em>.moment-2, the third as <em>label</em>.moment-3, etc.");
62 841 : }
63 :
64 5 : Moments::Moments( const vesselbase::VesselOptions& da) :
65 5 : StoreDataVessel(da)
66 : {
67 5 : ActionWithValue* a=dynamic_cast<ActionWithValue*>( getAction() );
68 5 : plumed_massert(a,"cannot create passable values as base action does not inherit from ActionWithValue");
69 :
70 5 : std::vector<std::string> moments=Tools::getWords(getAllInput(),"\t\n ,");
71 5 : Tools::interpretRanges(moments); unsigned nn;
72 11 : for(unsigned i=0; i<moments.size(); ++i) {
73 6 : a->addComponentWithDerivatives( "moment-" + moments[i] );
74 6 : a->componentIsNotPeriodic( "moment-" + moments[i] );
75 6 : value_out.push_back( a->copyOutput( a->getNumberOfComponents()-1 ) );
76 6 : Tools::convert( moments[i], nn );
77 6 : if( nn<2 ) error("moments are only possible for m>=2" );
78 6 : powers.push_back( nn ); std::string num; Tools::convert(powers[i],num);
79 11 : }
80 5 : }
81 :
82 229 : void Moments::resize() {
83 229 : StoreDataVessel::resize();
84 229 : }
85 :
86 5 : std::string Moments::description() {
87 10 : std::string descri, num;
88 5 : Tools::convert(powers[0],num);
89 5 : descri = "value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
90 6 : for(unsigned i=1; i<powers.size(); ++i) {
91 1 : Tools::convert(powers[i],num);
92 1 : descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
93 : }
94 10 : return descri;
95 : }
96 :
97 610 : void Moments::finish( const std::vector<double>& buffer ) {
98 610 : StoreDataVessel::finish( buffer );
99 :
100 610 : const double pi=3.141592653589793238462643383279502884197169399375105820974944592307;
101 610 : unsigned nvals=getAction()->getFullNumberOfTasks();
102 :
103 610 : double mean=0; Value myvalue;
104 610 : if( getAction()->isPeriodic() ) {
105 0 : std::string str_min, str_max; getAction()->retrieveDomain( str_min, str_max );
106 0 : double pfactor, min, max; Tools::convert(str_min,min); Tools::convert(str_max,max);
107 0 : pfactor = 2*pi / ( max-min ); myvalue.setDomain( str_min, str_max );
108 0 : double sinsum=0, cossum=0, val;
109 0 : for(unsigned i=0; i<nvals; ++i) { val=pfactor*( buffer[bufstart + i*nspace*vecsize+nspace] - min ); sinsum+=sin(val); cossum+=cos(val); }
110 0 : mean = 0.5 + atan2( sinsum / static_cast<double>( nvals ), cossum / static_cast<double>( nvals ) ) / (2*pi);
111 0 : mean = min + (max-min)*mean;
112 : } else {
113 610 : for(unsigned i=0; i<nvals; ++i) mean+=buffer[bufstart + i*nspace*vecsize+nspace];
114 610 : mean/=static_cast<double>( nvals ); myvalue.setNotPeriodic();
115 : }
116 :
117 1420 : for(unsigned npow=0; npow<powers.size(); ++npow) {
118 810 : double dev1=0;
119 810 : if( value_out[0]->getNumberOfDerivatives()>0 ) {
120 410 : for(unsigned i=0; i<nvals; ++i) dev1+=pow( myvalue.difference( mean, buffer[bufstart + i*nspace*vecsize+nspace] ), powers[npow] - 1 );
121 410 : dev1/=static_cast<double>( nvals );
122 : }
123 :
124 810 : double moment=0;
125 810 : MultiValue myvals( getNumberOfComponents(), getAction()->getNumberOfDerivatives() ); myvals.clearAll();
126 7710 : for(unsigned i=0; i<nvals; ++i) {
127 6900 : double tmp=myvalue.difference( mean, buffer[bufstart + i*nspace*vecsize+nspace] );
128 6900 : moment+=pow( tmp, powers[npow] );
129 6900 : if( value_out[npow]->getNumberOfDerivatives() ) {
130 4100 : double pref=pow( tmp, powers[npow] - 1 ) - dev1;
131 4100 : retrieveDerivatives( i, false, myvals );
132 164000 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
133 159900 : unsigned jatom=myvals.getActiveIndex(j);
134 159900 : value_out[npow]->addDerivative(jatom, pref*myvals.getDerivative( 1, jatom ) );
135 : }
136 4100 : myvals.clearAll();
137 : }
138 : }
139 810 : if( value_out[npow]->getNumberOfDerivatives()>0 ) value_out[npow]->chainRule( powers[npow] / static_cast<double>( nvals ) );
140 810 : value_out[npow]->set( moment / static_cast<double>( nvals ) );
141 1420 : }
142 610 : }
143 :
144 220 : bool Moments::applyForce( std::vector<double>& forces ) {
145 220 : std::vector<double> tmpforce( forces.size() );
146 220 : forces.assign(forces.size(),0.0); bool wasforced=false;
147 640 : for(unsigned i=0; i<value_out.size(); ++i) {
148 420 : if( value_out[i]->applyForce( tmpforce ) ) {
149 0 : wasforced=true;
150 0 : for(unsigned j=0; j<forces.size(); ++j) forces[j]+=tmpforce[j];
151 : }
152 : }
153 220 : return wasforced;
154 : }
155 :
156 : }
157 2523 : }
|
