Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-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 "AverageOnGrid.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace gridtools {
      26             : 
      27          65 : void AverageOnGrid::registerKeywords( Keywords& keys ) {
      28          65 :   HistogramOnGrid::registerKeywords( keys );
      29          65 : }
      30             : 
      31           8 : AverageOnGrid::AverageOnGrid( const vesselbase::VesselOptions& da ):
      32           8 :   HistogramOnGrid(da)
      33             : {
      34           8 :   arg_names.push_back( "density" );
      35           8 :   if( !discrete ) {
      36          22 :     for(unsigned i=0; i<dimension; ++i) arg_names.push_back( "ddensity_" + arg_names[i] );
      37           8 :     nper += (dimension+1);
      38             :   } else {
      39           0 :     nper += 1;
      40             :   }
      41           8 : }
      42             : 
      43    11825111 : void AverageOnGrid::accumulate( const unsigned& ipoint, const double& weight, const double& dens, const std::vector<double>& der, std::vector<double>& buffer ) const {
      44    11825111 :   buffer[bufstart+nper*ipoint] += weight*dens; buffer[ bufstart+nper*(ipoint+1) - (dimension+1) ] += dens;
      45    11825111 :   if( der.size()>0 ) {
      46    47300036 :     for(unsigned j=0; j<dimension; ++j) buffer[ bufstart+nper*ipoint + 1 + j ] += weight*der[j];
      47    47300036 :     for(unsigned j=0; j<dimension; ++j) buffer[ bufstart+nper*(ipoint+1) - dimension + j ] += der[j];
      48             :   }
      49    11825111 : }
      50             : 
      51      361702 : double AverageOnGrid::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
      52      361702 :   if( noAverage() ) return getDataElement( nper*ipoint + jelement);
      53             : 
      54      360502 :   if( jelement>=(nper-(dimension+1)) ) return getDataElement( nper*ipoint + jelement );
      55             : 
      56      360002 :   if( noderiv ) return getDataElement( nper*ipoint+jelement ) / getDataElement( nper*(1+ipoint) - 1);
      57             : 
      58             :   double rdenom = 1.0;
      59      360002 :   if( std::fabs(getDataElement( nper*(ipoint+1) -(dimension+1) ))>epsilon ) rdenom = 1. / getDataElement( nper*(ipoint+1) - (dimension+1) );
      60             : 
      61      360002 :   unsigned jderiv = jelement%(1+dimension);
      62      360002 :   if( jderiv==0 ) return rdenom*getDataElement( nper*ipoint+jelement );
      63             : 
      64      203082 :   unsigned jfloor = std::floor( jelement / (1+dimension) );
      65      203082 :   return rdenom*getDataElement( nper*ipoint+jelement ) - rdenom*rdenom*getDataElement(nper*ipoint+jfloor)*getDataElement(nper*(ipoint+1) - (dimension+1) + jderiv);
      66             : }
      67             : 
      68             : }
      69             : }