Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-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 "KernelFunctions.h"
      23             : #include "IFile.h"
      24             : #include <iostream>
      25             : #include <cmath>
      26             : 
      27             : namespace PLMD {
      28             : 
      29             : //+PLUMEDOC INTERNAL kernelfunctions
      30             : /*
      31             : Functions that are used to construct histograms
      32             : 
      33             : Constructing histograms is something you learned to do relatively early in life. You perform an experiment a number of times,
      34             : count the number of times each result comes up and then draw a bar graph that describes how often each of the results came up.
      35             : This only works when there are a finite number of possible results.  If the result a number between 0 and 1 the bar chart is
      36             : less easy to draw as there are as many possible results as there are numbers between zero and one - an infinite number.
      37             : To resolve this problem we replace probability, \f$P\f$ with probability density, \f$\pi\f$, and write the probability of getting
      38             : a number between \f$a\f$ and \f$b\f$ as:
      39             : 
      40             : \f[
      41             : P = \int_{a}^b \textrm{d}x \pi(x)
      42             : \f]
      43             : 
      44             : To calculate probability densities from a set of results we use a process called kernel density estimation.
      45             : Histograms are accumulated by adding up kernel functions, \f$K\f$, with finite spatial extent, that integrate to one.
      46             : These functions are centered on each of the \f$n\f$-dimensional data points, \f$\mathbf{x}_i\f$. The overall effect of this
      47             : is that each result we obtain in our experiments contributes to the probability density in a finite sized region of the space.
      48             : 
      49             : Expressing all this mathematically in kernel density estimation we write the probability density as:
      50             : 
      51             : \f[
      52             : \pi(\mathbf{x}) =  \sum_i K\left[ (\mathbf{x} - \mathbf{x}_i)^T \Sigma (\mathbf{x} - \mathbf{x}_i) \right]
      53             : \f]
      54             : 
      55             : where \f$\Sigma\f$ is an \f$n \times n\f$ matrix called the bandwidth that controls the spatial extent of
      56             : the kernel. Whenever we accumulate a histogram (e.g. in \ref HISTOGRAM or in \ref METAD) we use this
      57             : technique.
      58             : 
      59             : There is thus some flexibility in the particular function we use for \f$K[\mathbf{r}]\f$ in the above.
      60             : The following variants are available.
      61             : 
      62             : <table align=center frame=void width=95%% cellpadding=5%%>
      63             : <tr>
      64             : <td> TYPE </td> <td> FUNCTION </td>
      65             : </tr> <tr>
      66             : <td> gaussian </td> <td> \f$f(r) = \frac{1}{(2 \pi)^{n} \sqrt{|\Sigma^{-1}|}} \exp\left(-0.5 r^2 \right)\f$ </td>
      67             : </tr> <tr>
      68             : <td> truncated-gaussian </td> <td> \f$f(r) = \frac{1}{(2 \pi)^{n} \sqrt{|\Sigma^{-1}|} \left(\frac{\mathrm{erf}(-6.25/sqrt{2}) - \mathrm{erf}(-6.25/sqrt{2})}{2}\right)^n} \exp\left(-0.5 r^2 \right)\f$ </td>
      69             : </tr> <tr>
      70             : <td> triangular </td> <td> \f$f(r) = \frac{3}{V} ( 1 - | r | )H(1-|r|) \f$ </td>
      71             : </tr> <tr>
      72             : <td> uniform </td> <td> \f$f(r) = \frac{1}{V}H(1-|r|)\f$ </td>
      73             : </tr>
      74             : </table>
      75             : 
      76             : In the above \f$H(y)\f$ is a function that is equal to one when \f$y>0\f$ and zero when \f$y \le 0\f$. \f$n\f$ is
      77             : the dimensionality of the vector \f$\mathbf{x}\f$ and \f$V\f$ is the volume of an ellipse in an \f$n\f$ dimensional
      78             : space which is given by:
      79             : 
      80             : \f{eqnarray*}{
      81             : V &=& | \Sigma^{-1} | \frac{ \pi^{\frac{n}{2}} }{\left( \frac{n}{2} \right)! } \qquad \textrm{for even} \quad n \\
      82             : V &=& | \Sigma^{-1} | \frac{ 2^{\frac{n+1}{2}} \pi^{\frac{n-1}{2}} }{ n!! }
      83             : \f}
      84             : 
      85             : In \ref METAD the normalization constants are ignored so that the value of the function at \f$r=0\f$ is equal
      86             : to one.  In addition in \ref METAD we must be able to differentiate the bias in order to get forces.  This limits
      87             : the kernels we can use in this method.  Notice also that Gaussian kernels should have infinite support.  When used
      88             : with grids, however, they are assumed to only be non-zero over a finite range.  The difference between the
      89             : truncated-gaussian and regular gaussian is that the truncated gaussian is scaled so that its integral over the grid
      90             : is equal to one when it is normalized.  The integral of a regular gaussian when it is evaluated on a grid will be
      91             : slightly less that one because of the truncation of a function that should have infinite support.
      92             : */
      93             : //+ENDPLUMEDOC
      94             : 
      95           0 : KernelFunctions::KernelFunctions( const std::string& input ) {
      96             : 
      97           0 :   if(!dp2cutoffNoStretch()) {
      98           0 :     stretchA=dp2cutoffA;
      99           0 :     stretchB=dp2cutoffB;
     100             :   }
     101             : 
     102           0 :   std::vector<std::string> data=Tools::getWords(input);
     103           0 :   std::string name=data[0];
     104             :   data.erase(data.begin());
     105             : 
     106             :   std::vector<double> at;
     107           0 :   bool foundc = Tools::parseVector(data,"CENTER",at);
     108           0 :   if(!foundc) plumed_merror("failed to find center keyword in definition of kernel");
     109             :   std::vector<double> sig;
     110           0 :   bool founds = Tools::parseVector(data,"SIGMA",sig);
     111           0 :   if(!founds) plumed_merror("failed to find sigma keyword in definition of kernel");
     112             : 
     113           0 :   bool multi=false; Tools::parseFlag(data,"MULTIVARIATE",multi);
     114           0 :   bool vonmises=false; Tools::parseFlag(data,"VON-MISSES",vonmises);
     115           0 :   if( center.size()==1 && multi ) plumed_merror("one dimensional kernel cannot be multivariate");
     116           0 :   if( center.size()==1 && vonmises ) plumed_merror("one dimensional kernal cannot be von-misses");
     117           0 :   if( center.size()==1 && sig.size()!=1 ) plumed_merror("size mismatch between center size and sigma size");
     118           0 :   if( multi && center.size()>1 && sig.size()!=0.5*center.size()*(center.size()-1) ) plumed_merror("size mismatch between center size and sigma size");
     119           0 :   if( !multi && center.size()>1 && sig.size()!=center.size() ) plumed_merror("size mismatch between center size and sigma size");
     120             : 
     121             :   double h;
     122           0 :   bool foundh = Tools::parse(data,"HEIGHT",h);
     123           0 :   if( !foundh) h=1.0;
     124             : 
     125           0 :   if( multi ) setData( at, sig, name, "MULTIVARIATE", h );
     126           0 :   else if( vonmises ) setData( at, sig, name, "VON-MISSES", h );
     127           0 :   else setData( at, sig, name, "DIAGONAL", h );
     128           0 : }
     129             : 
     130        5563 : KernelFunctions::KernelFunctions( const std::vector<double>& at, const std::vector<double>& sig, const std::string& type, const std::string& mtype, const double& w ) {
     131             : 
     132        5563 :   if(!dp2cutoffNoStretch()) {
     133        5563 :     stretchA=dp2cutoffA;
     134        5563 :     stretchB=dp2cutoffB;
     135             :   }
     136             : 
     137        5563 :   setData( at, sig, type, mtype, w );
     138        5563 : }
     139             : 
     140           0 : KernelFunctions::KernelFunctions( const KernelFunctions* in ):
     141           0 :   dtype(in->dtype),
     142           0 :   ktype(in->ktype),
     143           0 :   center(in->center),
     144           0 :   width(in->width),
     145           0 :   height(in->height)
     146             : {
     147           0 :   if(!dp2cutoffNoStretch()) {
     148           0 :     stretchA=dp2cutoffA;
     149           0 :     stretchB=dp2cutoffB;
     150             :   }
     151           0 : }
     152             : 
     153        5563 : void KernelFunctions::setData( const std::vector<double>& at, const std::vector<double>& sig, const std::string& type, const std::string& mtype, const double& w ) {
     154             : 
     155        5563 :   height=w;
     156       16688 :   center.resize( at.size() ); for(unsigned i=0; i<at.size(); ++i) center[i]=at[i];
     157       21056 :   width.resize( sig.size() ); for(unsigned i=0; i<sig.size(); ++i) width[i]=sig[i];
     158        5563 :   if( mtype=="MULTIVARIATE" ) dtype=multi;
     159        1195 :   else if( mtype=="VON-MISSES" ) dtype=vonmises;
     160        1195 :   else if( mtype=="DIAGONAL" ) dtype=diagonal;
     161           0 :   else plumed_merror(mtype + " is not a valid metric type");
     162             : 
     163             :   // Setup the kernel type
     164       11126 :   if(type=="GAUSSIAN" || type=="gaussian" ) {
     165           0 :     ktype=gaussian;
     166       11126 :   } else if(type=="STRETCHED-GAUSSIAN" || type=="stretched-gaussian" ) {
     167        5563 :     ktype=stretchedgaussian;
     168           0 :   } else if(type=="TRUNCATED-GAUSSIAN" || type=="truncated-gaussian" ) {
     169           0 :     ktype=truncatedgaussian;
     170           0 :   } else if(type=="UNIFORM" || type=="uniform") {
     171           0 :     ktype=uniform;
     172           0 :   } else if(type=="TRIANGULAR" || type=="triangular") {
     173           0 :     ktype=triangular;
     174             :   } else {
     175           0 :     plumed_merror(type+" is an invalid kernel type\n");
     176             :   }
     177        5563 : }
     178             : 
     179           0 : void KernelFunctions::normalize( const std::vector<Value*>& myvals ) {
     180             : 
     181             :   double det=1.;
     182             :   unsigned ncv=ndim();
     183           0 :   if(dtype==diagonal) {
     184           0 :     for(unsigned i=0; i<width.size(); ++i) det*=width[i]*width[i];
     185           0 :   } else if(dtype==multi) {
     186           0 :     Matrix<double> mymatrix( getMatrix() ), myinv( ncv, ncv );
     187           0 :     Invert(mymatrix,myinv); double logd;
     188           0 :     logdet( myinv, logd );
     189           0 :     det=std::exp(logd);
     190             :   }
     191           0 :   if( dtype==diagonal || dtype==multi ) {
     192             :     double volume;
     193           0 :     if( ktype==gaussian || ktype==stretchedgaussian ) {
     194           0 :       volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 );
     195           0 :     } else if( ktype==truncatedgaussian ) {
     196             :       // This makes it so the gaussian integrates to one over the range over which it has support
     197             :       const double DP2CUTOFF=std::sqrt(6.25);
     198           0 :       volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 ) * pow( 0.5 * ( erf(DP2CUTOFF) - erf(-DP2CUTOFF) ), ncv);
     199           0 :     } else if( ktype==uniform || ktype==triangular ) {
     200           0 :       if( ncv%2==1 ) {
     201             :         double dfact=1;
     202           0 :         for(unsigned i=1; i<ncv; i+=2) dfact*=static_cast<double>(i);
     203           0 :         volume=( pow( pi, (ncv-1)/2 ) ) * ( pow( 2., (ncv+1)/2 ) ) / dfact;
     204             :       } else {
     205             :         double fact=1.;
     206           0 :         for(unsigned i=1; i<ncv/2; ++i) fact*=static_cast<double>(i);
     207           0 :         volume=pow( pi,ncv/2 ) / fact;
     208             :       }
     209           0 :       if(ktype==uniform) volume*=det;
     210           0 :       else if(ktype==triangular) volume*=det / 3.;
     211             :     } else {
     212           0 :       plumed_merror("not a valid kernel type");
     213             :     }
     214           0 :     height /= volume;
     215           0 :     return;
     216             :   }
     217           0 :   plumed_assert( dtype==vonmises && ktype==gaussian );
     218             :   // Now calculate determinant for aperiodic variables
     219             :   unsigned naper=0;
     220           0 :   for(unsigned i=0; i<ncv; ++i) {
     221           0 :     if( !myvals[i]->isPeriodic() ) naper++;
     222             :   }
     223             :   // Now construct sub matrix
     224             :   double volume=1;
     225           0 :   if( naper>0 ) {
     226             :     unsigned isub=0;
     227           0 :     Matrix<double> mymatrix( getMatrix() ), mysub( naper, naper );
     228           0 :     for(unsigned i=0; i<ncv; ++i) {
     229           0 :       if( myvals[i]->isPeriodic() ) continue;
     230             :       unsigned jsub=0;
     231           0 :       for(unsigned j=0; j<ncv; ++j) {
     232           0 :         if( myvals[j]->isPeriodic() ) continue;
     233           0 :         mysub( isub, jsub ) = mymatrix( i, j ); jsub++;
     234             :       }
     235           0 :       isub++;
     236             :     }
     237             :     Matrix<double> myisub( naper, naper ); double logd;
     238           0 :     Invert( mysub, myisub ); logdet( myisub, logd );
     239           0 :     det=std::exp(logd);
     240           0 :     volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 );
     241             :   }
     242             : 
     243             :   // Calculate volume of periodic variables
     244             :   unsigned nper=0;
     245           0 :   for(unsigned i=0; i<ncv; ++i) {
     246           0 :     if( myvals[i]->isPeriodic() ) nper++;
     247             :   }
     248             : 
     249             :   // Now construct sub matrix
     250           0 :   if( nper>0 ) {
     251             :     unsigned isub=0;
     252           0 :     Matrix<double> mymatrix( getMatrix() ),  mysub( nper, nper );
     253           0 :     for(unsigned i=0; i<ncv; ++i) {
     254           0 :       if( !myvals[i]->isPeriodic() ) continue;
     255             :       unsigned jsub=0;
     256           0 :       for(unsigned j=0; j<ncv; ++j) {
     257           0 :         if( !myvals[j]->isPeriodic() ) continue;
     258           0 :         mysub( isub, jsub ) = mymatrix( i, j ); jsub++;
     259             :       }
     260           0 :       isub++;
     261             :     }
     262             :     Matrix<double>  eigvec( nper, nper );
     263           0 :     std::vector<double> eigval( nper );
     264           0 :     diagMat( mysub, eigval, eigvec );
     265             :     unsigned iper=0; volume=1;
     266           0 :     for(unsigned i=0; i<ncv; ++i) {
     267           0 :       if( myvals[i]->isPeriodic() ) {
     268           0 :         volume *= myvals[i]->getMaxMinusMin()*Tools::bessel0(eigval[iper])*std::exp(-eigval[iper]);
     269           0 :         iper++;
     270             :       }
     271             :     }
     272             :   }
     273           0 :   height /= volume;
     274             : }
     275             : 
     276       10033 : double KernelFunctions::getCutoff( const double& width ) const {
     277             :   const double DP2CUTOFF=6.25;
     278       10033 :   if( ktype==gaussian || ktype==truncatedgaussian || ktype==stretchedgaussian ) return std::sqrt(2.0*DP2CUTOFF)*width;
     279           0 :   else if(ktype==triangular ) return width;
     280           0 :   else if(ktype==uniform) return width;
     281           0 :   else plumed_merror("No valid kernel type");
     282             :   return 0.0;
     283             : }
     284             : 
     285        5017 : std::vector<double> KernelFunctions::getContinuousSupport( ) const {
     286             :   unsigned ncv=ndim();
     287        5017 :   std::vector<double> support( ncv );
     288        5017 :   if(dtype==diagonal) {
     289        1946 :     for(unsigned i=0; i<ncv; ++i) support[i]=getCutoff(width[i]);
     290        4368 :   } else if(dtype==multi) {
     291        4368 :     Matrix<double> mymatrix( getMatrix() ), myinv( ncv,ncv );
     292        4368 :     Invert(mymatrix,myinv);
     293        4368 :     Matrix<double> myautovec(ncv,ncv); std::vector<double> myautoval(ncv);
     294        4368 :     diagMat(myinv,myautoval,myautovec);
     295             :     double maxautoval=0.;
     296             :     unsigned ind_maxautoval=0;
     297       13104 :     for (unsigned i=0; i<ncv; i++) {
     298        8736 :       if(myautoval[i]>maxautoval) {maxautoval=myautoval[i]; ind_maxautoval=i;}
     299             :     }
     300       13104 :     for(unsigned i=0; i<ncv; ++i) {
     301        8736 :       double extent=std::fabs(std::sqrt(maxautoval)*myautovec(i,ind_maxautoval));
     302        8736 :       support[i]=getCutoff( extent );
     303             :     }
     304             :   } else {
     305           0 :     plumed_merror("cannot find support if metric is not multi or diagonal type");
     306             :   }
     307        5017 :   return support;
     308             : }
     309             : 
     310        3375 : std::vector<unsigned> KernelFunctions::getSupport( const std::vector<double>& dx ) const {
     311        3375 :   plumed_assert( ndim()==dx.size() );
     312        3375 :   std::vector<unsigned> support( dx.size() );
     313        3375 :   std::vector<double> vv=getContinuousSupport( );
     314       10124 :   for(unsigned i=0; i<dx.size(); ++i) support[i]=static_cast<unsigned>(ceil( vv[i]/dx[i] ));
     315        3375 :   return support;
     316             : }
     317             : 
     318     1024305 : double KernelFunctions::evaluate( const std::vector<Value*>& pos, std::vector<double>& derivatives, bool usederiv, bool doInt, double lowI_, double uppI_) const {
     319             :   plumed_dbg_assert( pos.size()==ndim() && derivatives.size()==ndim() );
     320             : #ifndef NDEBUG
     321             :   if( usederiv ) plumed_massert( ktype!=uniform, "step function can not be differentiated" );
     322             : #endif
     323     1024305 :   if(doInt) {
     324             :     plumed_dbg_assert(center.size()==1);
     325           0 :     if(pos[0]->get()<lowI_) pos[0]->set(lowI_);
     326           0 :     if(pos[0]->get()>uppI_) pos[0]->set(uppI_);
     327             :   }
     328             :   double r2=0;
     329     1024305 :   if(dtype==diagonal) {
     330      161311 :     for(unsigned i=0; i<ndim(); ++i) {
     331      107510 :       derivatives[i]=-pos[i]->difference( center[i] ) / width[i];
     332      107510 :       r2+=derivatives[i]*derivatives[i];
     333      107510 :       derivatives[i] /= width[i];
     334             :     }
     335      970504 :   } else if(dtype==multi) {
     336      970504 :     Matrix<double> mymatrix( getMatrix() );
     337     2911512 :     for(unsigned i=0; i<mymatrix.nrows(); ++i) {
     338     1941008 :       double dp_i, dp_j; derivatives[i]=0;
     339     1941008 :       dp_i=-pos[i]->difference( center[i] );
     340     5823024 :       for(unsigned j=0; j<mymatrix.ncols(); ++j) {
     341     3882016 :         if(i==j) dp_j=dp_i;
     342     1941008 :         else dp_j=-pos[j]->difference( center[j] );
     343             : 
     344     3882016 :         derivatives[i]+=mymatrix(i,j)*dp_j;
     345     3882016 :         r2+=dp_i*dp_j*mymatrix(i,j);
     346             :       }
     347             :     }
     348           0 :   } else if(dtype==vonmises) {
     349           0 :     std::vector<double> costmp( ndim() ), sintmp( ndim() ), sinout( ndim(), 0.0 );
     350           0 :     for(unsigned i=0; i<ndim(); ++i) {
     351           0 :       if( pos[i]->isPeriodic() ) {
     352           0 :         sintmp[i]=std::sin( 2.*pi*(pos[i]->get() - center[i])/pos[i]->getMaxMinusMin() );
     353           0 :         costmp[i]=std::cos( 2.*pi*(pos[i]->get() - center[i])/pos[i]->getMaxMinusMin() );
     354             :       } else {
     355           0 :         sintmp[i]=pos[i]->get() - center[i];
     356           0 :         costmp[i]=1.0;
     357             :       }
     358             :     }
     359             : 
     360           0 :     Matrix<double> mymatrix( getMatrix() );
     361           0 :     for(unsigned i=0; i<mymatrix.nrows(); ++i) {
     362           0 :       derivatives[i]=0;
     363           0 :       if( pos[i]->isPeriodic() ) {
     364           0 :         r2+=2*( 1 - costmp[i] )*mymatrix(i,i);
     365             :       } else {
     366           0 :         r2+=sintmp[i]*sintmp[i]*mymatrix(i,i);
     367             :       }
     368           0 :       for(unsigned j=0; j<mymatrix.ncols(); ++j) {
     369           0 :         if( i!=j ) sinout[i]+=mymatrix(i,j)*sintmp[j];
     370             :       }
     371           0 :       derivatives[i] = mymatrix(i,i)*sintmp[i] + sinout[i]*costmp[i];
     372           0 :       if( pos[i]->isPeriodic() ) derivatives[i] *= (2*pi/pos[i]->getMaxMinusMin());
     373             :     }
     374           0 :     for(unsigned i=0; i<sinout.size(); ++i) r2+=sintmp[i]*sinout[i];
     375             :   }
     376             :   double kderiv, kval;
     377     1024305 :   if(ktype==gaussian || ktype==truncatedgaussian) {
     378           0 :     kval=height*std::exp(-0.5*r2); kderiv=-kval;
     379     1024305 :   } else if(ktype==stretchedgaussian) {
     380     1024305 :     auto dp=0.5*r2;
     381     1024305 :     if(dp<dp2cutoff) {
     382      630478 :       auto ee=std::exp(-0.5*r2);
     383      630478 :       kval=height*(stretchA*ee+stretchB);
     384      630478 :       kderiv=-height*stretchA*ee;
     385             :     } else {
     386             :       kval=0.0;
     387             :       kderiv=0.0;
     388             :     }
     389             :   } else {
     390           0 :     double r=std::sqrt(r2);
     391           0 :     if(ktype==triangular) {
     392           0 :       if( r<1.0 ) {
     393             :         if(r==0) kderiv=0;
     394           0 :         kderiv=-1; kval=height*( 1. - std::fabs(r) );
     395             :       } else {
     396             :         kval=0.; kderiv=0.;
     397             :       }
     398           0 :     } else if(ktype==uniform) {
     399             :       kderiv=0.;
     400           0 :       if(r<1.0) kval=height;
     401             :       else kval=0;
     402             :     } else {
     403           0 :       plumed_merror("Not a valid kernel type");
     404             :     }
     405           0 :     kderiv*=height / r ;
     406             :   }
     407     3072823 :   for(unsigned i=0; i<ndim(); ++i) derivatives[i]*=kderiv;
     408     1024305 :   if(doInt) {
     409           0 :     if((pos[0]->get() <= lowI_ || pos[0]->get() >= uppI_) && usederiv ) for(unsigned i=0; i<ndim(); ++i)derivatives[i]=0;
     410             :   }
     411     1024305 :   return kval;
     412             : }
     413             : 
     414        4471 : std::unique_ptr<KernelFunctions> KernelFunctions::read( IFile* ifile, const bool& cholesky, const std::vector<std::string>& valnames ) {
     415             :   double h;
     416        8942 :   if( !ifile->scanField("height",h) ) return NULL;;
     417             : 
     418        4471 :   std::string sss; ifile->scanField("multivariate",sss);
     419       12321 :   std::string ktype="stretched-gaussian"; if( ifile->FieldExist("kerneltype") ) ifile->scanField("kerneltype",ktype);
     420        8839 :   plumed_massert( sss=="false" || sss=="true" || sss=="von-misses", "multivariate flag must be either false, true or von-misses");
     421             : 
     422             :   // Read the position of the center
     423        4471 :   std::vector<double> cc( valnames.size() );
     424       13412 :   for(unsigned i=0; i<valnames.size(); ++i) ifile->scanField(valnames[i],cc[i]);
     425             : 
     426             :   std::vector<double> sig;
     427        4471 :   if( sss=="false" ) {
     428         103 :     sig.resize( valnames.size() );
     429         308 :     for(unsigned i=0; i<valnames.size(); ++i) {
     430         205 :       ifile->scanField("sigma_"+valnames[i],sig[i]);
     431         205 :       if( !cholesky ) sig[i]=std::sqrt(sig[i]);
     432             :     }
     433             :     return Tools::make_unique<KernelFunctions>(cc, sig, ktype, "DIAGONAL", h);
     434             :   }
     435             : 
     436        4368 :   unsigned ncv=valnames.size();
     437        4368 :   sig.resize( (ncv*(ncv+1))/2 );
     438             :   Matrix<double> upper(ncv,ncv), lower(ncv,ncv), mymult( ncv, ncv ), invmatrix(ncv,ncv);
     439       13104 :   for(unsigned i=0; i<ncv; ++i) {
     440       21840 :     for(unsigned j=0; j<ncv-i; j++) {
     441       26208 :       ifile->scanField("sigma_" +valnames[j+i] + "_" + valnames[j], lower(j+i,j) );
     442       13104 :       upper(j,j+i)=lower(j+i,j); mymult(j+i,j)=mymult(j,j+i)=lower(j+i,j);
     443             :     }
     444             :   }
     445        4368 :   if( cholesky ) mult(lower,upper,mymult);
     446        4368 :   Invert( mymult, invmatrix );
     447             :   unsigned k=0;
     448       13104 :   for(unsigned i=0; i<ncv; i++) {
     449       21840 :     for(unsigned j=i; j<ncv; j++) { sig[k]=invmatrix(i,j); k++; }
     450             :   }
     451        4368 :   if( sss=="true" ) return Tools::make_unique<KernelFunctions>(cc, sig, ktype, "MULTIVARIATE", h);
     452             :   return Tools::make_unique<KernelFunctions>( cc, sig, ktype, "VON-MISSES", h );
     453             : }
     454             : 
     455             : }