diff --git a/src/ves/LinearBasisSetExpansion.cpp b/src/ves/LinearBasisSetExpansion.cpp
index e17e14852fede47c3743ecc4d1e40841c4526c3f..343072de9fcf69209db73cbb1f19e2f58229048d 100644
--- a/src/ves/LinearBasisSetExpansion.cpp
+++ b/src/ves/LinearBasisSetExpansion.cpp
@@ -468,6 +468,31 @@ void LinearBasisSetExpansion::getBasisSetValues(const std::vector<double>& args_
 }
 
 
+double LinearBasisSetExpansion::getBasisSetValue(const std::vector<double>& args_values, const size_t index, std::vector<BasisFunctions*>& basisf_pntrs_in, CoeffsVector* coeffs_pntr_in) {
+  unsigned int nargs = args_values.size();
+  plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
+  plumed_assert(basisf_pntrs_in.size()==nargs);
+
+  std::vector<double> args_values_trsfrm(nargs);
+  std::vector< std::vector <double> > bf_values;
+  //
+  for(unsigned int k=0; k<nargs; k++) {
+    std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
+    std::vector<double> tmp_der(tmp_val.size());
+    bool inside=true;
+    basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
+    bf_values.push_back(tmp_val);
+  }
+  //
+  std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(index);
+  double bf_value=1.0;
+  for(unsigned int k=0; k<nargs; k++) {
+    bf_value*=bf_values[k][indices[k]];
+  }
+  return bf_value;
+}
+
+
 void LinearBasisSetExpansion::setupUniformTargetDistribution() {
   std::vector< std::vector <double> > bf_integrals(0);
   std::vector<double> targetdist_averages(ncoeffs_,0.0);
diff --git a/src/ves/LinearBasisSetExpansion.h b/src/ves/LinearBasisSetExpansion.h
index a7872d216c371f3fae82fafe0a54dc0c71cf959c..e4e581b1a60b46da4b51ea3e64a3f1ab87466bf1 100644
--- a/src/ves/LinearBasisSetExpansion.h
+++ b/src/ves/LinearBasisSetExpansion.h
@@ -125,6 +125,10 @@ public:
   //
   static void getBasisSetValues(const std::vector<double>&, std::vector<double>&, std::vector<BasisFunctions*>&, CoeffsVector*, Communicator* comm_in=NULL);
   void getBasisSetValues(const std::vector<double>&, std::vector<double>&, const bool parallel=true);
+  //
+  static double getBasisSetValue(const std::vector<double>&, const size_t, std::vector<BasisFunctions*>&, CoeffsVector*);
+  double getBasisSetValue(const std::vector<double>&, const size_t);
+  double getBasisSetConstant();
   // Bias grid and output stuff
   void setupBiasGrid(const bool usederiv=false);
   void updateBiasGrid();
@@ -224,6 +228,19 @@ void LinearBasisSetExpansion::getBasisSetValues(const std::vector<double>& args_
 }
 
 
+inline
+double LinearBasisSetExpansion::getBasisSetValue(const std::vector<double>& args_values, const size_t basisset_index) {
+  return getBasisSetValue(args_values,basisset_index,basisf_pntrs_, bias_coeffs_pntr_);
+}
+
+
+inline
+double LinearBasisSetExpansion::getBasisSetConstant() {
+  std::vector<double> args_dummy(nargs_,0.0);
+  return getBasisSetValue(args_dummy,0,basisf_pntrs_, bias_coeffs_pntr_);
+}
+
+
 }
 
 }