Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "SphericalKDEGridTools.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace gridtools {
      26             : 
      27          46 : void SphericalKDEGridTools::registerKeywords( Keywords& keys ) {
      28          46 :   keys.add("compulsory","CONCENTRATION","the concentration parameter for the Von Mises-Fisher distributions");
      29          46 :   keys.add("compulsory","HEIGHTS","1.0","this keyword takes the label of an action that calculates a vector of values. The elements of this vector "
      30             :            "are used as weights for the Gaussians.");
      31          46 :   keys.add("compulsory","GRID_BIN","the number of points on the fibonacci sphere at which the density should be evaluated");
      32          46 : }
      33             : 
      34          13 : void SphericalKDEGridTools::readBandwidthAndHeight( const UniversalVonMisses& params, ActionWithArguments* action ) {
      35             :   // Read in the concentration parameters
      36             :   std::string von_misses_concentration;
      37          13 :   action->parse("CONCENTRATION",von_misses_concentration);
      38          13 :   KDEHelper<VonMissesKernelParams,UniversalVonMisses,SphericalKDEGridTools>::readKernelParameters( von_misses_concentration, action, "_vmconcentration", true );
      39          13 :   action->log.printf("  getting concentration parameters from %s \n", von_misses_concentration.c_str() );
      40             :   // Read in the heights
      41             :   std::string weight_str;
      42          13 :   action->parse("HEIGHTS",weight_str);
      43          26 :   KDEHelper<VonMissesKernelParams,UniversalVonMisses,SphericalKDEGridTools>::readKernelParameters( weight_str, action, "_volumes", false );
      44          13 :   if( (action->getPntrToArgument(0))->getNumberOfValues()==1 ) {
      45           0 :     action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO ARG=" + weight_str + "," + action->getLabel() + "_vmconcentration FUNC=x*y/(4*pi*sinh(y))" ), false );
      46             :   } else {
      47          52 :     action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO ARG=" + weight_str + "," + action->getLabel() + "_vmconcentration FUNC=x*y/(4*pi*sinh(y)) MASK=" + weight_str ), false );
      48             :   }
      49          13 :   KDEHelper<VonMissesKernelParams,UniversalVonMisses,SphericalKDEGridTools>::addArgument( action->getLabel() + "_heights", action );
      50          13 :   action->log.printf("  getting heights from %s \n", weight_str.c_str() );
      51          13 : }
      52             : 
      53          13 : void SphericalKDEGridTools::readGridParameters( SphericalKDEGridTools& g, ActionWithArguments* action, GridCoordinatesObject& gridobject, std::vector<std::size_t>& shape ) {
      54          13 :   if( shape.size()!=3 ) {
      55           0 :     action->error("should have three coordinates in input to this action");
      56             :   }
      57          13 :   action->parse("GRID_BIN",g.nbins);
      58          13 :   action->log.printf("  setting number of bins to %zu \n", g.nbins );
      59          13 :   std::vector<bool> ipbc( 3, false );
      60          14 :   gridobject.setup( "fibonacci", ipbc, g.nbins, 0 );
      61          12 :   shape[0]=g.nbins;
      62          12 :   shape[1]=shape[2]=1;
      63          12 : }
      64             : 
      65          11 : void SphericalKDEGridTools::getDiscreteSupport( const SphericalKDEGridTools& g, const UniversalVonMisses& p, const View<const double>& shape, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject ) {
      66          11 :   plumed_assert( nneigh.size()==gridobject.getDimension() );
      67          11 :   std::vector<bool> ipbc( 3, false );
      68          11 :   double fib_cutoff = std::log( epsilon / (shape[0]/(4*pi*sinh(shape[0]))) ) / shape[0];   // The shape here is the concentration of the fisher kernel
      69          22 :   gridobject.setup( "fibonacci", ipbc, gridobject.getNumberOfPoints(), fib_cutoff );
      70          11 : }
      71             : 
      72      132180 : void SphericalKDEGridTools::getNeighbors( const UniversalVonMisses& p, View<double> at, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors ) {
      73      132180 :   gridobject.getNeighbors( View<const double>(at.data(),at.size()), nneigh, num_neighbors, neighbors );
      74      132180 : }
      75             : 
      76             : }
      77             : }