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           8 :   arg_names.push_back( "density" );
      34           8 :   if( !discrete ) {
      35          22 :     for(unsigned i=0; i<dimension; ++i) {
      36          28 :       arg_names.push_back( "ddensity_" + arg_names[i] );
      37             :     }
      38           8 :     nper += (dimension+1);
      39             :   } else {
      40           0 :     nper += 1;
      41             :   }
      42           8 : }
      43             : 
      44    11825111 : void AverageOnGrid::accumulate( const unsigned& ipoint, const double& weight, const double& dens, const std::vector<double>& der, std::vector<double>& buffer ) const {
      45    11825111 :   buffer[bufstart+nper*ipoint] += weight*dens;
      46    11825111 :   buffer[ bufstart+nper*(ipoint+1) - (dimension+1) ] += dens;
      47    11825111 :   if( der.size()>0 ) {
      48    47300036 :     for(unsigned j=0; j<dimension; ++j) {
      49    35474925 :       buffer[ bufstart+nper*ipoint + 1 + j ] += weight*der[j];
      50             :     }
      51    47300036 :     for(unsigned j=0; j<dimension; ++j) {
      52    35474925 :       buffer[ bufstart+nper*(ipoint+1) - dimension + j ] += der[j];
      53             :     }
      54             :   }
      55    11825111 : }
      56             : 
      57      361734 : double AverageOnGrid::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
      58      361734 :   if( noAverage() ) {
      59        1200 :     return getDataElement( nper*ipoint + jelement);
      60             :   }
      61             : 
      62      360534 :   if( jelement>=(nper-(dimension+1)) ) {
      63         500 :     return getDataElement( nper*ipoint + jelement );
      64             :   }
      65             : 
      66      360034 :   if( noderiv ) {
      67           0 :     return getDataElement( nper*ipoint+jelement ) / getDataElement( nper*(1+ipoint) - 1);
      68             :   }
      69             : 
      70             :   double rdenom = 1.0;
      71      360034 :   if( std::fabs(getDataElement( nper*(ipoint+1) -(dimension+1) ))>epsilon ) {
      72      317549 :     rdenom = 1. / getDataElement( nper*(ipoint+1) - (dimension+1) );
      73             :   }
      74             : 
      75      360034 :   unsigned jderiv = jelement%(1+dimension);
      76      360034 :   if( jderiv==0 ) {
      77      156928 :     return rdenom*getDataElement( nper*ipoint+jelement );
      78             :   }
      79             : 
      80      203106 :   unsigned jfloor = std::floor( jelement / (1+dimension) );
      81      203106 :   return rdenom*getDataElement( nper*ipoint+jelement ) - rdenom*rdenom*getDataElement(nper*ipoint+jfloor)*getDataElement(nper*(ipoint+1) - (dimension+1) + jderiv);
      82             : }
      83             : 
      84             : }
      85             : }