Working pamm collective variables and tests based on real data

regtest/multicolvar/rt-pamm-periodic/2D-testc-0.75.pammp

+#! FIELDS weight center-phi center-psi width-phi-phi width-phi-psi width-psi-phi width-psi-psi      
+      2.97197455E-0001     -1.91983118E+0000      2.25029540E+0000      2.45960237E-0001     -1.30615381E-0001     -1.30615381E-0001      2.40239117E-0001  
+      2.29131448E-0002      1.39809354E+0000      9.54585380E-0002      9.61755708E-0002     -3.55657919E-0002     -3.55657919E-0002      1.06147253E-0001  
+      5.06676398E-0001     -1.09648066E+0000     -7.17867907E-0001      1.40523052E-0001     -1.05385552E-0001     -1.05385552E-0001      1.63290557E-0001  
+      8.38235835E-0003      1.34746362E+0000     -2.92912242E+0000      1.50601708E-0001     -1.16056708E-0001     -1.16056708E-0001      2.58187042E-0001  
+      1.34976025E-0001     -1.12488647E+0000      2.46856186E+0000      6.47837871E-0002     -5.48993853E-0002     -5.48993853E-0002      1.28240303E-0001  
+      2.62404899E-0002      9.93277439E-0001      6.63990207E-0001      3.77594608E-0002     -2.13908228E-0002     -2.13908228E-0002      7.02294006E-0002  
+      3.61412967E-0003      9.92586372E-0001     -2.28743740E+0000      6.13474371E-0002      9.32283874E-0003      9.32283874E-0003      4.90062255E-0002  

regtest/multicolvar/rt-pamm-periodic/colvar.reference

+#! FIELDS time p.mean
+ 0.000000   0.9911
+ 1.000000   0.9980
+ 2.000000   0.9981
+ 3.000000   0.9959
+ 4.000000   0.9981
+ 5.000000   0.9954

regtest/multicolvar/rt-pamm-periodic/config

+type=driver
+# this is to test a different name
+arg="--plumed plumed.dat --ixyz M1d_open.xyz --box 55.598,55.569,55.569 --length-units A"

regtest/multicolvar/rt-pamm-periodic/plumed.dat

+MOLINFO MOLTYPE=protein STRUCTURE=M1d.pdb
+
+TORSIONS ...
+   ATOMS1=@psi-2
+   ATOMS2=@psi-3
+   ATOMS3=@psi-4
+   ATOMS4=@psi-5
+   ATOMS5=@psi-6
+   ATOMS6=@psi-7
+   ATOMS7=@psi-8
+   ATOMS8=@psi-9
+   ATOMS9=@psi-10
+   ATOMS10=@psi-11
+   ATOMS11=@psi-12
+   ATOMS12=@psi-13
+   ATOMS13=@psi-14
+   ATOMS14=@psi-15
+   ATOMS15=@psi-16
+   ATOMS16=@psi-17
+   ATOMS17=@psi-18
+   ATOMS18=@psi-19
+   ATOMS19=@psi-20
+   ATOMS20=@psi-21
+   ATOMS21=@psi-22
+   LABEL=psi
+... TORSIONS
+
+TORSIONS ...
+   ATOMS1=@phi-2
+   ATOMS2=@phi-3
+   ATOMS3=@phi-4
+   ATOMS4=@phi-5
+   ATOMS5=@phi-6
+   ATOMS6=@phi-7
+   ATOMS7=@phi-8
+   ATOMS8=@phi-9
+   ATOMS9=@phi-10
+   ATOMS10=@phi-11
+   ATOMS11=@phi-12
+   ATOMS12=@phi-13
+   ATOMS13=@phi-14
+   ATOMS14=@phi-15
+   ATOMS15=@phi-16
+   ATOMS16=@phi-17
+   ATOMS17=@phi-18
+   ATOMS18=@phi-19
+   ATOMS19=@phi-20
+   ATOMS20=@phi-21
+   ATOMS21=@phi-22
+   LABEL=phi
+... TORSIONS
+
+DUMPMULTICOLVAR DATA=psi FILE=psi.xyz
+DUMPMULTICOLVAR DATA=phi FILE=phi.xyz
+
+PAMM DATA=phi,psi CLUSTERS=2D-testc-0.75.pammp MEAN LABEL=p
+# PAMM DATA=phi,psi CLUSTERS=2D-testc-0.75.pammp MEAN NUMERICAL_DERIVATIVES LABEL=pn
+
+PRINT ARG=p.* FILE=colvar FMT=%8.4f
+DUMPDERIVATIVES ARG=p.* FILE=deriv FMT=%8.4f
+

regtest/multicolvar/rt-pamm/clusters.dat

-#! FIELDS weight center-d1 center-d2 width-d1-d1 width-d1-d2 width-d2-d1 width-d2-d2
- 0.4 1.2 1.2 10.0 0.0 0.0 10.0
- 0.6 1.5 1.5 10.0 0.0 0.0 10.0

regtest/multicolvar/rt-pamm/colvar.reference

-#! FIELDS time p.mean
- 0.000000   0.5794
- 0.050000   0.5255
- 0.100000   0.4830
- 0.150000   0.4394
- 0.200000   0.4340

regtest/multicolvar/rt-pamm/config

-type=driver
-# this is to test a different name
-arg="--plumed plumed.dat --trajectory-stride 10 --timestep 0.005 --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt=%10.6f --pdb test.pdb"
-extra_files="../../trajectories/trajectory.xyz"

regtest/multicolvar/rt-pamm/deriv.reference

-#! FIELDS time parameter p.mean pn.mean
- 0.000000 0   0.1847   0.1847
- 0.000000 1  -0.0024  -0.0024
- 0.000000 2   0.1629   0.1629
- 0.000000 3  -0.1847  -0.1847
- 0.000000 4   0.0024   0.0024
- 0.000000 5  -0.1629  -0.1629
- 0.000000 6   0.1831   0.1831
- 0.000000 7  -0.0100  -0.0100
- 0.000000 8   0.1546   0.1546
- 0.000000 9  -0.1831  -0.1831
- 0.000000 10   0.0100   0.0100
- 0.000000 11  -0.1546  -0.1546
- 0.000000 12   0.1721   0.1721
- 0.000000 13  -0.0007  -0.0007
- 0.000000 14   0.1732   0.1732
- 0.000000 15  -0.1721  -0.1721
- 0.000000 16   0.0007   0.0007
- 0.000000 17  -0.1732  -0.1732
- 0.000000 18  -0.1745  -0.1745
- 0.000000 19  -0.0056  -0.0056
- 0.000000 20   0.1737   0.1737
- 0.000000 21   0.1745   0.1745
- 0.000000 22   0.0056   0.0056
- 0.000000 23  -0.1737  -0.1737
- 0.000000 24   0.1698   0.1698
- 0.000000 25  -0.0005  -0.0005
- 0.000000 26   0.1694   0.1694
- 0.000000 27  -0.1698  -0.1698
- 0.000000 28   0.0005   0.0005
- 0.000000 29  -0.1694  -0.1694
- 0.000000 30  -0.1772  -0.1772
- 0.000000 31   0.0023   0.0023
- 0.000000 32   0.1679   0.1679
- 0.000000 33   0.1772   0.1772
- 0.000000 34  -0.0023  -0.0023
- 0.000000 35  -0.1679  -0.1679
- 0.000000 36   0.9491   0.0000
- 0.000000 37  -0.0097   0.0000
- 0.000000 38   0.2914   0.0000
- 0.000000 39  -0.0097   0.0000
- 0.000000 40   0.0007   0.0000
- 0.000000 41  -0.0136   0.0000
- 0.000000 42   0.2914   0.0000
- 0.000000 43  -0.0136   0.0000
- 0.000000 44   0.8454   0.0000
- 0.050000 0   0.1943   0.1943
- 0.050000 1  -0.0015  -0.0015
- 0.050000 2   0.1568   0.1568
- 0.050000 3  -0.1943  -0.1943
- 0.050000 4   0.0015   0.0015
- 0.050000 5  -0.1568  -0.1568
- 0.050000 6   0.1973   0.1973
- 0.050000 7  -0.0166  -0.0166
- 0.050000 8   0.1458   0.1458
- 0.050000 9  -0.1973  -0.1973
- 0.050000 10   0.0166   0.0166
- 0.050000 11  -0.1458  -0.1458
- 0.050000 12   0.1729   0.1729
- 0.050000 13  -0.0016  -0.0016
- 0.050000 14   0.1795   0.1795
- 0.050000 15  -0.1729  -0.1729
- 0.050000 16   0.0016   0.0016
- 0.050000 17  -0.1795  -0.1795
- 0.050000 18  -0.1795  -0.1795
- 0.050000 19  -0.0116  -0.0116
- 0.050000 20   0.1732   0.1732
- 0.050000 21   0.1795   0.1795
- 0.050000 22   0.0116   0.0116
- 0.050000 23  -0.1732  -0.1732
- 0.050000 24   0.1758   0.1758
- 0.050000 25  -0.0051  -0.0051
- 0.050000 26   0.1717   0.1717
- 0.050000 27  -0.1758  -0.1758
- 0.050000 28   0.0051   0.0051
- 0.050000 29  -0.1717  -0.1717
- 0.050000 30  -0.1823  -0.1823
- 0.050000 31   0.0038   0.0038
- 0.050000 32   0.1698   0.1698
- 0.050000 33   0.1823   0.1823
- 0.050000 34  -0.0038  -0.0038
- 0.050000 35  -0.1698  -0.1698
- 0.050000 36   1.0351   0.0000
- 0.050000 37  -0.0169   0.0000
- 0.050000 38   0.2981   0.0000
- 0.050000 39  -0.0169   0.0000
- 0.050000 40   0.0023   0.0000
- 0.050000 41  -0.0262   0.0000
- 0.050000 42   0.2981   0.0000
- 0.050000 43  -0.0262   0.0000
- 0.050000 44   0.8476   0.0000
- 0.100000 0   0.1893   0.1893
- 0.100000 1  -0.0015  -0.0015
- 0.100000 2   0.1459   0.1459
- 0.100000 3  -0.1893  -0.1893
- 0.100000 4   0.0015   0.0015
- 0.100000 5  -0.1459  -0.1459
- 0.100000 6   0.2035   0.2035
- 0.100000 7  -0.0183  -0.0183
- 0.100000 8   0.1405   0.1405
- 0.100000 9  -0.2035  -0.2035
- 0.100000 10   0.0183   0.0183
- 0.100000 11  -0.1405  -0.1405
- 0.100000 12   0.1705   0.1705
- 0.100000 13  -0.0075  -0.0075
- 0.100000 14   0.1826   0.1826
- 0.100000 15  -0.1705  -0.1705
- 0.100000 16   0.0075   0.0075
- 0.100000 17  -0.1826  -0.1826
- 0.100000 18  -0.1778  -0.1778
- 0.100000 19  -0.0121  -0.0121
- 0.100000 20   0.1592   0.1592
- 0.100000 21   0.1778   0.1778
- 0.100000 22   0.0121   0.0121
- 0.100000 23  -0.1592  -0.1592
- 0.100000 24   0.1807   0.1807
- 0.100000 25  -0.0035  -0.0035
- 0.100000 26   0.1698   0.1698
- 0.100000 27  -0.1807  -0.1807
- 0.100000 28   0.0035   0.0035
- 0.100000 29  -0.1698  -0.1698
- 0.100000 30  -0.1791  -0.1791
- 0.100000 31   0.0058   0.0058
- 0.100000 32   0.1742   0.1742
- 0.100000 33   0.1791   0.1791
- 0.100000 34  -0.0058  -0.0058
- 0.100000 35  -0.1742  -0.1742
- 0.100000 36   1.0702   0.0000
- 0.100000 37  -0.0240   0.0000
- 0.100000 38   0.3042   0.0000
- 0.100000 39  -0.0240   0.0000
- 0.100000 40   0.0031   0.0000
- 0.100000 41  -0.0299   0.0000
- 0.100000 42   0.3042   0.0000
- 0.100000 43  -0.0299   0.0000
- 0.100000 44   0.8328   0.0000
- 0.150000 0   0.1673   0.1673
- 0.150000 1   0.0007   0.0007
- 0.150000 2   0.1254   0.1254
- 0.150000 3  -0.1673  -0.1673
- 0.150000 4  -0.0007  -0.0007
- 0.150000 5  -0.1254  -0.1254
- 0.150000 6   0.2028   0.2028
- 0.150000 7  -0.0222  -0.0222
- 0.150000 8   0.1415   0.1415
- 0.150000 9  -0.2028  -0.2028
- 0.150000 10   0.0222   0.0222
- 0.150000 11  -0.1415  -0.1415
- 0.150000 12   0.1698   0.1698
- 0.150000 13  -0.0135  -0.0135
- 0.150000 14   0.1824   0.1824
- 0.150000 15  -0.1698  -0.1698
- 0.150000 16   0.0135   0.0135
- 0.150000 17  -0.1824  -0.1824
- 0.150000 18  -0.1585  -0.1585
- 0.150000 19  -0.0068  -0.0068
- 0.150000 20   0.1361   0.1361
- 0.150000 21   0.1585   0.1585
- 0.150000 22   0.0068   0.0068
- 0.150000 23  -0.1361  -0.1361
- 0.150000 24   0.1838   0.1838
- 0.150000 25  -0.0032  -0.0032
- 0.150000 26   0.1669   0.1669
- 0.150000 27  -0.1838  -0.1838
- 0.150000 28   0.0032   0.0032
- 0.150000 29  -0.1669  -0.1669
- 0.150000 30  -0.1779  -0.1779
- 0.150000 31   0.0108   0.0108
- 0.150000 32   0.1747   0.1747
- 0.150000 33   0.1779   0.1779
- 0.150000 34  -0.0108  -0.0108
- 0.150000 35  -0.1747  -0.1747
- 0.150000 36   1.0570   0.0000
- 0.150000 37  -0.0405   0.0000
- 0.150000 38   0.3045   0.0000
- 0.150000 39  -0.0405   0.0000
- 0.150000 40   0.0045   0.0000
- 0.150000 41  -0.0275   0.0000
- 0.150000 42   0.3045   0.0000
- 0.150000 43  -0.0275   0.0000
- 0.150000 44   0.8045   0.0000
- 0.200000 0   0.1561   0.1561
- 0.200000 1   0.0004   0.0004
- 0.200000 2   0.1109   0.1109
- 0.200000 3  -0.1561  -0.1561
- 0.200000 4  -0.0004  -0.0004
- 0.200000 5  -0.1109  -0.1109
- 0.200000 6   0.1962   0.1962
- 0.200000 7  -0.0274  -0.0274
- 0.200000 8   0.1460   0.1460
- 0.200000 9  -0.1962  -0.1962
- 0.200000 10   0.0274   0.0274
- 0.200000 11  -0.1460  -0.1460
- 0.200000 12   0.1730   0.1730
- 0.200000 13  -0.0085  -0.0085
- 0.200000 14   0.1798   0.1798
- 0.200000 15  -0.1730  -0.1730
- 0.200000 16   0.0085   0.0085
- 0.200000 17  -0.1798  -0.1798
- 0.200000 18  -0.1438  -0.1438
- 0.200000 19  -0.0036  -0.0036
- 0.200000 20   0.1264   0.1264
- 0.200000 21   0.1438   0.1438
- 0.200000 22   0.0036   0.0036
- 0.200000 23  -0.1264  -0.1264
- 0.200000 24   0.1800   0.1800
- 0.200000 25  -0.0002  -0.0002
- 0.200000 26   0.1678   0.1678
- 0.200000 27  -0.1800  -0.1800
- 0.200000 28   0.0002   0.0002
- 0.200000 29  -0.1678  -0.1678
- 0.200000 30  -0.1797  -0.1797
- 0.200000 31   0.0132   0.0132
- 0.200000 32   0.1728   0.1728
- 0.200000 33   0.1797   0.1797
- 0.200000 34  -0.0132  -0.0132
- 0.200000 35  -0.1728  -0.1728
- 0.200000 36   1.0251   0.0000
- 0.200000 37  -0.0442   0.0000
- 0.200000 38   0.2964   0.0000
- 0.200000 39  -0.0442   0.0000
- 0.200000 40   0.0054   0.0000
- 0.200000 41  -0.0205   0.0000
- 0.200000 42   0.2964   0.0000
- 0.200000 43  -0.0205   0.0000
- 0.200000 44   0.7851   0.0000

regtest/multicolvar/rt-pamm/plumed.dat

-DISTANCES ATOMS1=1,2 ATOMS2=3,40 ATOMS3=5,6 LABEL=d1
-DISTANCES ATOMS1=7,8 ATOMS2=9,10 ATOMS3=11,12 LABEL=d2
-
-PAMM DATA=d1,d2 CLUSTERS=clusters.dat MEAN LABEL=p
-PAMM DATA=d1,d2 CLUSTERS=clusters.dat MEAN NUMERICAL_DERIVATIVES LABEL=pn
-
-PRINT ARG=p.* FILE=colvar FMT=%8.4f
-DUMPDERIVATIVES ARG=p.*,pn.* FILE=deriv FMT=%8.4f 
-

regtest/multicolvar/rt-pamm/test.pdb

-ATOM      1  Ar                  0.000   0.000   0.000  1.00 +0.50
-ATOM      2  Ar                  0.000   0.000   0.000  1.00 +0.50
-ATOM      3  Ar                  0.000   0.000   0.000  1.00 +0.50
-ATOM      4  Ar                  0.000   0.000   0.000  1.00 +0.50
-ATOM      5  Ar                  0.000   0.000   0.000  1.00 +0.50
-ATOM      6  Ar                  0.000   0.000   0.000  1.00 -0.50
-ATOM      7  Ar                  0.000   0.000   0.000  1.00 -0.50
-ATOM      8  Ar                  0.000   0.000   0.000  1.00 -0.50
-ATOM      9  Ar                  0.000   0.000   0.000  1.00 -0.50
-ATOM     10  Ar                  0.000   0.000   0.000  1.00 -0.50
-END

src/multicolvar/PAMM.cpp

@@ -91,23 +91,26 @@ MultiColvarFunction(ao)
    }
 
    std::string filename; parse("CLUSTERS",filename);
-   IFile ifile; 
+   unsigned ncolv = getNumberOfBaseMultiColvars();
+   IFile ifile; Matrix<double> covar( ncolv, ncolv ), icovar( ncolv, ncolv );
    if( !ifile.FileExist(filename) ) error("could not find file named " + filename);
    ifile.open(filename); ifile.allowIgnoredFields(); double weight, wij;
-   unsigned ncolv = getNumberOfBaseMultiColvars(); 
    std::vector<double> center( ncolv );
    std::vector<double> sig( 0.5*ncolv*(ncolv-1) + ncolv );   
    while( ifile.scanField("weight",weight) ) {
        for(unsigned i=0;i<center.size();++i){
           ifile.scanField("center-"+getBaseMultiColvar(i)->getLabel(),center[i]);
        }   
-       unsigned k=0;
        for(unsigned i=0;i<center.size();++i){
            for(unsigned j=0;j<center.size();++j){
-              ifile.scanField("width-"+getBaseMultiColvar(i)->getLabel()+"-" + getBaseMultiColvar(j)->getLabel(),wij );
-              if( i<=j) sig[k++]=wij;
+              ifile.scanField("width-"+getBaseMultiColvar(i)->getLabel()+"-" + getBaseMultiColvar(j)->getLabel(),covar(i,j) );
            }
        }  
+       Invert( covar, icovar );
+       unsigned k=0; 
+       for(unsigned i=0;i<ncolv;i++){
+           for(unsigned j=i;j<ncolv;j++){ sig[k]=icovar(i,j); k++; }
+       }
        ifile.scanField(); 
        // std::string ktype; ifile.scanField("kerneltype",ktype);
        kernels.push_back( KernelFunctions( center, sig, "GAUSSIAN", "VON-MISES", weight ) );

src/tools/KernelFunctions.cpp

@@ -151,12 +151,38 @@ void KernelFunctions::normalize( const std::vector<Value*>& myvals ){
      Invert(mymatrix,myinv); double logd;
      logdet( myinv, logd );
      det=std::exp(logd);
-  } else if(dtype==vonmises){
-     unsigned naper=0; 
-     for(unsigned i=0;i<ncv;++i){
-         if( !myvals[i]->isPeriodic() ) naper++;
+  } 
+  if( dtype==diagonal || dtype==multi ){
+     double volume;
+     if( ktype==gaussian ){
+        volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 );
+     } else if( ktype==uniform || ktype==triangular ){
+        if( ncv%2==1 ){
+           double dfact=1;
+           for(unsigned i=1;i<ncv;i+=2) dfact*=static_cast<double>(i);
+           volume=( pow( pi, (ncv-1)/2 ) ) * ( pow( 2., (ncv+1)/2 ) ) / dfact;
+        } else {
+           double fact=1.;
+           for(unsigned i=1;i<ncv/2;++i) fact*=static_cast<double>(i);
+           volume=pow( pi,ncv/2 ) / fact;
+        }
+        if(ktype==uniform) volume*=det;
+        else if(ktype==triangular) volume*=det / 3.;
+     } else {
+        plumed_merror("not a valid kernel type");
      }
-     // Now construct sub matrix
+     height /= volume;
+     return;
+  }
+  plumed_assert( dtype==vonmises && ktype==gaussian );
+  // Now calculate determinant for aperiodic variables
+  unsigned naper=0; 
+  for(unsigned i=0;i<ncv;++i){
+      if( !myvals[i]->isPeriodic() ) naper++;
+  }
+  // Now construct sub matrix
+  double volume=1;
+  if( naper>0 ){
      unsigned isub=0;
      Matrix<double> mymatrix( getMatrix() ), mysub( naper, naper );
      for(unsigned i=0;i<ncv;++i){
@@ -171,26 +197,40 @@ void KernelFunctions::normalize( const std::vector<Value*>& myvals ){
      Matrix<double> myisub( naper, naper ); double logd;
      Invert( mysub, myisub ); logdet( myisub, logd );
      det=std::exp(logd);
-  }
-  double volume;
-  if( ktype==gaussian ){
      volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 );
-  } else if( ktype==uniform || ktype==triangular ){
-     if( ncv%2==1 ){
-        double dfact=1;
-        for(unsigned i=1;i<ncv;i+=2) dfact*=static_cast<double>(i);
-        volume=( pow( pi, (ncv-1)/2 ) ) * ( pow( 2., (ncv+1)/2 ) ) / dfact;
-     } else {
-        double fact=1.;
-        for(unsigned i=1;i<ncv/2;++i) fact*=static_cast<double>(i);
-        volume=pow( pi,ncv/2 ) / fact;
+  }
+
+  // Calculate volume of periodic variables 
+  unsigned nper=0;
+  for(unsigned i=0;i<ncv;++i){
+     if( myvals[i]->isPeriodic() ) nper++;
+  }
+
+  // Now construct sub matrix
+  if( nper>0 ){
+     unsigned isub=0;
+     Matrix<double> mymatrix( getMatrix() ),  mysub( nper, nper ); 
+     for(unsigned i=0;i<ncv;++i){
+         if( !myvals[i]->isPeriodic() ) continue;
+         unsigned jsub=0;
+         for(unsigned j=0;j<ncv;++j){
+             if( !myvals[j]->isPeriodic() ) continue;
+             mysub( isub, jsub ) = mymatrix( i, j ); jsub++;
+         }
+         isub++;
      }
-     if(ktype==uniform) volume*=det;
-     else if(ktype==triangular) volume*=det / 3.;
-  } else {
-     plumed_merror("not a valid kernel type");
-  } 
-  height /= volume;  
+     Matrix<double>  eigvec( nper, nper ); 
+     std::vector<double> eigval( nper ); 
+     diagMat( mysub, eigval, eigvec ); 
+     unsigned iper=0; volume=1;
+     for(unsigned i=0;i<ncv;++i){
+        if( myvals[i]->isPeriodic() ){ 
+            volume *= myvals[i]->getMaxMinusMin()*Tools::bessel0(eigval[iper])*std::exp(-eigval[iper]);
+            iper++;
+        }
+     }
+  }
+  height /= volume;
 }
 
 double KernelFunctions::getCutoff( const double& width ) const {

src/tools/MolDataClass.cpp

@@ -67,6 +67,9 @@ bool MolDataClass::allowedResidue( const std::string& type, const std::string& r
       else if(residuename=="HSE") return true; // HIS-E charmm
       else if(residuename=="HIP") return true; // HIS-P amber
       else if(residuename=="HSP") return true; // HIS-P charmm
+// Weird amino acids
+      else if(residuename=="NLE") return true;
+      else if(residuename=="SFO") return true;
       else return false; 
   } else if( type=="dna" ){
       if(residuename=="DA") return true;
@@ -202,7 +205,7 @@ void MolDataClass::specialSymbol( const std::string& type, const std::string& sy
         numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum+1,chainid));
         numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum+1,chainid));
       } else if( name=="chi1" && !isTerminalGroup("protein",resname) ){
-        if ( resname=="GLY" || resname=="ALA" ) plumed_merror("chi-1 is not defined for Alanine and Glycine");
+        if ( resname=="GLY" || resname=="ALA" || resname=="SFO" ) plumed_merror("chi-1 is not defined for Alanine, Glycine and SFO");
         numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
         numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
         numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CB",resnum,chainid));

src/tools/Tools.cpp

@@ -322,4 +322,15 @@ std::string Tools::extension(const std::string&s){
   return ext;
 }
 
+double Tools::bessel0( const double& val ){
+ if (fabs(val)<3.75){
+      double y = Tools::fastpow( val/3.75, 2 );
+      return 1 + y*(3.5156229 +y*(3.0899424 + y*(1.2067492+y*(0.2659732+y*(0.0360768+y*0.0045813)))));
+  }
+  double ax=fabs(val), y=3.75/ax, bx=std::exp(ax)/sqrt(ax);
+  ax=0.39894228+y*(0.01328592+y*(0.00225319+y*(-0.00157565+y*(0.00916281+y*(-0.02057706+y*(0.02635537+y*(-0.01647633+y*0.00392377)))))));
+  return ax*bx;
+}
+
+
 }

src/tools/Tools.h

@@ -121,6 +121,8 @@ public:
   static std::string extension(const std::string&);
 /// Fast int power
   static double fastpow(double base,int exp);
+/// Modified 0th-order Bessel function of the first kind 
+  static double bessel0(const double& val);
 };
 
 template <class T>