Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-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 "vesselbase/VesselRegister.h"
      23             : #include "vesselbase/FunctionVessel.h"
      24             : #include "vesselbase/ActionWithVessel.h"
      25             : #include "multicolvar/ActionVolume.h"
      26             : #include "VectorMultiColvar.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace crystallization {
      30             : 
      31             : class VectorMean : public vesselbase::FunctionVessel {
      32             : private:
      33             :   unsigned nder;
      34             : public:
      35             :   static void registerKeywords( Keywords& keys );
      36             :   static void reserveKeyword( Keywords& keys );
      37             :   explicit VectorMean( const vesselbase::VesselOptions& da );
      38             :   std::string value_descriptor();
      39             :   void resize();
      40             :   void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
      41             :   void finish( const std::vector<double>& buffer );
      42             : };
      43             : 
      44       13791 : PLUMED_REGISTER_VESSEL(VectorMean,"VMEAN")
      45             : 
      46           6 : void VectorMean::registerKeywords( Keywords& keys ) {
      47           6 :   vesselbase::FunctionVessel::registerKeywords(keys);
      48           6 : }
      49             : 
      50        4595 : void VectorMean::reserveKeyword( Keywords& keys ) {
      51        9190 :   keys.reserve("vessel","VMEAN","calculate the norm of the mean vector.");
      52        9190 :   keys.addOutputComponent("vmean","VMEAN","the norm of the mean vector. The output component can be referred to elsewhere in the input "
      53             :                           "file by using the label.vmean");
      54        4595 : }
      55             : 
      56           6 : VectorMean::VectorMean( const vesselbase::VesselOptions& da ) :
      57             :   FunctionVessel(da),
      58           6 :   nder(0) {
      59           6 : }
      60             : 
      61           6 : std::string VectorMean::value_descriptor() {
      62           6 :   return "the norm of the mean vector";
      63             : }
      64             : 
      65          23 : void VectorMean::resize() {
      66          23 :   unsigned ncomp=getAction()->getNumberOfQuantities() - 2;
      67             : 
      68          23 :   if( getAction()->derivativesAreRequired() ) {
      69          16 :     nder=getAction()->getNumberOfDerivatives();
      70          16 :     resizeBuffer( (1+nder)*(ncomp+1) );
      71          16 :     getFinalValue()->resizeDerivatives( nder );
      72             :   } else {
      73           7 :     nder=0;
      74           7 :     resizeBuffer(ncomp+1);
      75             :   }
      76          23 : }
      77             : 
      78       14712 : void VectorMean::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
      79       14712 :   unsigned ncomp=getAction()->getNumberOfQuantities()-2;
      80             : 
      81       14712 :   double weight=myvals.get(0);
      82             :   plumed_dbg_assert( weight>=getTolerance() );
      83       14712 :   buffer[bufstart] += weight;
      84       61240 :   for(unsigned i=0; i<ncomp; ++i) {
      85       46528 :     buffer[bufstart + (1+i)*(1+nder)] += weight*myvals.get(2+i);
      86             :   }
      87       14712 :   if( !getAction()->derivativesAreRequired() ) {
      88             :     return;
      89             :   }
      90             : 
      91       14712 :   if( diffweight ) {
      92          29 :     myvals.chainRule( 0, 0, 1, 0, 1.0, bufstart, buffer );
      93             :   }
      94       61240 :   for(unsigned i=0; i<ncomp; ++i) {
      95       46528 :     double colvar=myvals.get(2+i);
      96       46528 :     myvals.chainRule( 2+i, 1+i, 1, 0, weight, bufstart, buffer );
      97       46528 :     if( diffweight ) {
      98         754 :       myvals.chainRule( 0, 1+i, 1, 0, colvar, bufstart, buffer );
      99             :     }
     100             :   }
     101             :   return;
     102             : }
     103             : 
     104        1828 : void VectorMean::finish( const std::vector<double>& buffer ) {
     105        1828 :   unsigned ncomp=getAction()->getNumberOfQuantities()-2;
     106        1828 :   double sum=0, ww=buffer[bufstart];
     107        7358 :   for(unsigned i=0; i<ncomp; ++i) {
     108        5530 :     double tmp = buffer[bufstart+(nder+1)*(i+1)] / ww;
     109        5530 :     sum+=tmp*tmp;
     110             :   }
     111        1828 :   setOutputValue( sqrt(sum) );
     112        1828 :   if( !getAction()->derivativesAreRequired() ) {
     113             :     return;
     114             :   }
     115             : 
     116             :   Value* fval=getFinalValue();
     117        1828 :   double tw = 1.0 / sqrt(sum);
     118        7358 :   for(unsigned icomp=0; icomp<ncomp; ++icomp) {
     119        5530 :     double tmp = buffer[bufstart + (icomp+1)*(1+nder)] / ww;
     120        5530 :     unsigned bstart = bufstart + (1+icomp)*(nder+1) + 1;
     121      460636 :     for(unsigned jder=0; jder<nder; ++jder) {
     122      455106 :       fval->addDerivative( jder, (tw*tmp/ww)*( buffer[bstart + jder] - tmp*buffer[bufstart + 1 + jder] ) );
     123             :     }
     124             :   }
     125             : }
     126             : 
     127             : }
     128             : }