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       10427 : 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        3473 : void Moments::reserveKeyword( Keywords& keys ) {
      61        6946 :   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        6946 :   keys.reset_style("MOMENTS","vessel");
      68        6946 :   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        3473 : }
      72             : 
      73           8 : Moments::Moments( const vesselbase::VesselOptions& da) :
      74           8 :   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 = "moment-";
      81          24 :   parse("COMPONENT",mycomponent); parseVector("MOMENTS",moments);
      82          10 :   if( getNumericalLabel()!=0 ) { std::string numstr; Tools::convert( mycomponent, numstr); valstr = "moment-" + numstr + "-"; }
      83          14 :   if( moments.size()==0 ) moments=Tools::getWords(getAllInput(),"\t\n ,");
      84           8 :   Tools::interpretRanges(moments); unsigned nn;
      85          20 :   for(unsigned i=0; i<moments.size(); ++i) {
      86          24 :     a->addComponentWithDerivatives( valstr + moments[i] );
      87          24 :     a->componentIsNotPeriodic( valstr + moments[i] );
      88          12 :     value_out.push_back( a->copyOutput( a->getNumberOfComponents()-1 ) );
      89          12 :     Tools::convert( moments[i], nn );
      90          12 :     if( nn<2 ) error("moments are only possible for m>=2" );
      91          12 :     powers.push_back( nn ); std::string num; Tools::convert(powers[i],num);
      92             :   }
      93           8 : }
      94             : 
      95          28 : void Moments::resize() {
      96             :   resizeBuffer(0);
      97          28 :   if( getAction()->derivativesAreRequired() ) {
      98          56 :     for(unsigned i=0; i<value_out.size(); ++i) value_out[i]->resizeDerivatives( getAction()->getNumberOfDerivatives() );
      99             :   }
     100          28 : }
     101             : 
     102           8 : std::string Moments::description() {
     103             :   std::string descri, num;
     104           8 :   Tools::convert(powers[0],num);
     105           8 :   if( getNumericalLabel()==0 ) {
     106          12 :     descri = "value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
     107           8 :     for(unsigned i=1; i<powers.size(); ++i) {
     108           2 :       Tools::convert(powers[i],num);
     109           4 :       descri = descri + "\n  value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
     110             :     }
     111             :   } else {
     112           2 :     std::string numlab; Tools::convert( mycomponent, numlab );
     113           4 :     descri = "value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
     114           4 :     for(unsigned i=1; i<powers.size(); ++i) {
     115           2 :       Tools::convert(powers[i],num);
     116           4 :       descri = descri + "\n  value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
     117             :     }
     118             :   }
     119           8 :   return descri;
     120             : }
     121             : 
     122         822 : void Moments::finish( const std::vector<double>& buffer ) {
     123             :   const double pi=3.141592653589793238462643383279502884197169399375105820974944592307;
     124             :   unsigned nvals=getAction()->getFullNumberOfTasks();
     125         822 :   std::vector<double>  myvalues( getAction()->getNumberOfQuantities() );
     126             : 
     127         822 :   double mean=0; Value myvalue;
     128         822 :   if( getAction()->isPeriodic() ) {
     129           0 :     std::string str_min, str_max; getAction()->retrieveDomain( str_min, str_max );
     130           0 :     double pfactor, min, max; Tools::convert(str_min,min); Tools::convert(str_max,max);
     131           0 :     pfactor = 2*pi / ( max-min ); myvalue.setDomain( str_min, str_max );
     132             :     double sinsum=0, cossum=0, val;
     133           0 :     for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
     134           0 :       mystash->retrieveSequentialValue( i, false, myvalues );
     135           0 :       val=pfactor*( myvalues[mycomponent] - min );
     136           0 :       sinsum+=std::sin(val); cossum+=std::cos(val);
     137             :     }
     138           0 :     mean = 0.5 + std::atan2( sinsum / static_cast<double>( nvals ), cossum / static_cast<double>( nvals ) ) / (2*pi);
     139           0 :     mean = min + (max-min)*mean;
     140             :   } else {
     141        7758 :     for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
     142        6936 :       mystash->retrieveSequentialValue( i, false, myvalues ); mean+=myvalues[mycomponent];
     143             :     }
     144         822 :     mean/=static_cast<double>( nvals ); myvalue.setNotPeriodic();
     145             :   }
     146             : 
     147        2056 :   for(unsigned npow=0; npow<powers.size(); ++npow) {
     148             :     double dev1=0;
     149        1234 :     if( value_out[0]->getNumberOfDerivatives()>0 ) {
     150        4606 :       for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
     151        4172 :         mystash->retrieveSequentialValue( i, false, myvalues );
     152        4172 :         dev1+=pow( myvalue.difference( mean, myvalues[mycomponent] ), powers[npow] - 1 );
     153             :       }
     154         434 :       dev1/=static_cast<double>( nvals );
     155             :     }
     156             : 
     157             :     double moment=0;
     158        1234 :     MultiValue myvals( getAction()->getNumberOfQuantities(), getAction()->getNumberOfDerivatives() ); myvals.clearAll();
     159       11006 :     for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
     160        9772 :       mystash->retrieveSequentialValue( i, false, myvalues );
     161        9772 :       double tmp=myvalue.difference( mean, myvalues[mycomponent] );
     162        9772 :       moment+=pow( tmp, powers[npow] );
     163        9772 :       if( value_out[npow]->getNumberOfDerivatives() ) {
     164        4172 :         double pref=pow( tmp, powers[npow] - 1 ) - dev1;
     165        4172 :         mystash->retrieveDerivatives( i, false, myvals );
     166      166880 :         for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
     167             :           unsigned jatom=myvals.getActiveIndex(j);
     168      162708 :           value_out[npow]->addDerivative(jatom, pref*myvals.getDerivative( mycomponent, jatom ) );
     169             :         }
     170        4172 :         myvals.clearAll();
     171             :       }
     172             :     }
     173        1234 :     if( value_out[npow]->getNumberOfDerivatives()>0 ) value_out[npow]->chainRule( powers[npow] / static_cast<double>( nvals ) );
     174        1234 :     value_out[npow]->set( moment / static_cast<double>( nvals ) );
     175        1234 :   }
     176        1644 : }
     177             : 
     178         432 : bool Moments::applyForce( std::vector<double>& forces ) {
     179         432 :   std::vector<double> tmpforce( forces.size() );
     180         432 :   forces.assign(forces.size(),0.0); bool wasforced=false;
     181        1276 :   for(unsigned i=0; i<value_out.size(); ++i) {
     182         844 :     if( value_out[i]->applyForce( tmpforce ) ) {
     183             :       wasforced=true;
     184           0 :       for(unsigned j=0; j<forces.size(); ++j) forces[j]+=tmpforce[j];
     185             :     }
     186             :   }
     187         432 :   return wasforced;
     188             : }
     189             : 
     190             : }
     191             : }