Added command to calculate multicolvar averages etc in a spherical region centered on an atom

regtest/multicolvar/rt-cylinder/Makefile

+include ../../scripts/test.make

regtest/multicolvar/rt-cylinder/colvar.reference

+#! FIELDS time d2 d2s
+ 0.000000   0.1249   3.1925
+ 0.050000   0.0308   4.4595
+ 0.100000   0.1379   4.2236
+ 0.150000   0.0731   3.4469
+ 0.200000   0.2866   3.3730
+ 0.250000   0.2583   3.9249
+ 0.300000   0.1894   4.7553
+ 0.350000   0.4102   5.8226
+ 0.400000   0.2946   2.7438
+ 0.450000   0.2556   3.3862
+ 0.500000   0.0240   3.2564
+ 0.550000   0.1111   3.0024
+ 0.600000   0.1547   3.5411
+ 0.650000   0.0653   3.6699
+ 0.700000   0.2409   4.3339
+ 0.750000   0.2259   3.9310
+ 0.800000   0.0927   3.3248
+ 0.850000   0.2652   3.5731
+ 0.900000   0.1454   4.1791
+ 0.950000   0.2917   3.0881

regtest/multicolvar/rt-cylinder/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=%8.4f"

regtest/multicolvar/rt-cylinder/gentraj.cpp

+#include <stdlib.h>
+#include <math.h>
+#include <iostream>
+#include <vector>
+
+double randomn() { return rand()/static_cast<double>(RAND_MAX); }
+
+int main() {
+
+  std::vector<double> tot(8,0), dtot(8,0); bool dens=true;
+  unsigned nmols=200, dir; std::vector<double> com(3), dist(3), a1(3), a2(3), len(nmols);
+  for(unsigned i=0;i<20;++i){
+      if(dens) std::cout<<1+nmols<<std::endl;
+      else std::cout<<1+2*nmols<<std::endl;
+      std::cout<<10.0<<" "<<10.0<<" "<<10.0<<std::endl;
+      std::cout<<"Ar "<<5.0<<" "<<5.0<<" "<<5.0<<std::endl;
+      for(unsigned j=0;j<nmols;++j){
+         com[0]=10.0*randomn(); com[1]=10.0*randomn(); com[2]=10.0*randomn();
+         for(unsigned k=0;k<3;++k) dist[k]=fabs(com[k]-5.0);
+         len[j]=randomn(); dir=floor( 3*randomn() );
+         if( dist[0]<1.0 ){ tot[1]+=1.0; dtot[1]+=len[j]; }
+         if( dist[1]<1.0 ){ tot[2]+=1.0; dtot[2]+=len[j]; }
+         if( dist[2]<1.0 ){ tot[3]+=1.0; dtot[3]+=len[j]; }
+         if( dist[0]<1.0 && dist[1]<1.0 ){ tot[4]+=1.0; dtot[4]+=len[j]; }
+         if( dist[0]<1.0 && dist[2]<1.0 ){ tot[5]+=1.0; dtot[5]+=len[j]; }
+         if( dist[1]<1.0 && dist[2]<1.0 ){ tot[6]+=1.0; dtot[6]+=len[j]; }
+         if( dist[0]<1.0 && dist[1]<1.0 && dist[2]<1.0 ){ tot[7]+=1.0; dtot[7]+=len[j]; }
+         a1[0]=com[0]; a1[1]=com[1]; a1[2]=com[2];
+         a2[0]=com[0]; a2[1]=com[1]; a2[2]=com[2];
+         a1[dir]-=0.5*len[j]; a2[dir]+=0.5*len[j];
+         if(dens){
+            std::cout<<"Ar "<<com[0]<<" "<<com[1]<<" "<<com[2]<<std::endl;
+         } else {
+            std::cout<<"Ar "<<a1[0]<<" "<<a1[1]<<" "<<a1[2]<<std::endl;
+            std::cout<<"Ar "<<a2[0]<<" "<<a2[1]<<" "<<a2[2]<<std::endl;
+         }
+      }
+      for(unsigned i=1;i<8;++i){
+         if(dens) std::cerr<<"CV "<<i<<" "<<tot[i]<<std::endl; 
+         else std::cerr<<"CV "<<i<<dtot[i] / static_cast<double>( nmols )<<std::endl;
+         tot[i]=dtot[i]=0.0;
+      }
+  }
+
+  return 1;
+}

regtest/multicolvar/rt-cylinder/plumed.dat

+c1: COM ATOMS=1-200 
+
+DENSITY SPECIES=2-201 LABEL=d1
+
+d2: INCYLINDER ATOM=c1 DATA=d1 DIRECTION=Z RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-0.1 UPPER=0.1
+
+# DENSITY SPECIES=2-201 LABEL=d1n
+# d2n: INCYLINDER ATOM=c1 DATA=d1n DIRECTION=Z RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-0.1 UPPER=0.1 NUMERICAL_DERIVATIVES
+
+d2s: INSPHERE ATOM=c1 DATA=d1 RADIUS={TANH R_0=1.5}
+
+PRINT ARG=d2,d2s FILE=colvar FMT=%8.4f 
+DUMPDERIVATIVES ARG=d2,d2s FILE=derivatives FMT=%8.4f
+
+

src/multicolvar/VolumeInSphere.cpp

+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2015 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed-code.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "tools/Pbc.h"
+#include "tools/SwitchingFunction.h"
+#include "ActionVolume.h"
+
+//+PLUMEDOC MCOLVARF INSPHERE
+/*
+This quantity can be used to calculate functions of the distribution of collective 
+variables for the atoms that lie in a particular, user-specified part of of the cell.
+
+Each of the base quantities calculated by a multicolvar can can be assigned to a particular point in three 
+dimensional space. For example, if we have the coordination numbers for all the atoms in the
+system each coordination number can be assumed to lie on the position of the central atom. 
+Because each base quantity can be assigned to a particular point in space we can calculate functions of the
+distribution of base quantities in a particular part of the box by using:
+
+\f[
+\overline{s}_{\tau} = \frac{ \sum_i f(s_i) w(x_i,y_i,z_i) }{ \sum_i w(x_i,y_i,z_i) }  
+\f]  
+
+where the sum is over the collective variables, \f$s_i\f$, each of which can be thought to be at \f$ (x_i,y_i,z_i)\f$.
+The function \f$ w(x_i,y_i,z_i) \f$ measures whether or not the system is in the subregion of interest. It
+is equal to:
+
+\f[
+w(x_i,y_i,z_i) =  
+\f]
+
+where \f$\sigma\f$ is a \ref switchingfunction.
+The function \f$(s_i)\f$ can be any of the usual LESS_THAN, MORE_THAN, WITHIN etc that are used in all other multicolvars.
+
+When INCYLINDER is used with the \ref DENSITY action the number of atoms in the specified region is calculated  
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace multicolvar {
+
+class VolumeInSphere : public ActionVolume {
+private:
+  Vector origin;
+  SwitchingFunction switchingFunction;
+public:
+  static void registerKeywords( Keywords& keys );
+  explicit VolumeInSphere(const ActionOptions& ao);
+  void setupRegions();
+  double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const ;
+}; 
+
+PLUMED_REGISTER_ACTION(VolumeInSphere,"INSPHERE")
+
+void VolumeInSphere::registerKeywords( Keywords& keys ){
+  ActionVolume::registerKeywords( keys );
+  keys.add("atoms","ATOM","the atom whose vicinity we are interested in examining");
+  keys.add("compulsory","RADIUS","the switching function that tells us the extent of the sphereical region of interest");
+  keys.remove("SIGMA");
+}
+
+VolumeInSphere::VolumeInSphere(const ActionOptions& ao):
+Action(ao),
+ActionVolume(ao)
+{
+  std::vector<AtomNumber> atom; 
+  parseAtomList("ATOM",atom);
+  if( atom.size()!=1 ) error("should only be one atom specified");
+  log.printf("  center of sphere is at position of atom : %d\n",atom[0].serial() );
+
+  std::string sw, errors; parse("RADIUS",sw);
+  if(sw.length()==0) error("missing RADIUS keyword");
+  switchingFunction.set(sw,errors);
+  if( errors.length()!=0 ) error("problem reading RADIUS keyword : " + errors );
+  log.printf("  radius of sphere is given by %s \n", ( switchingFunction.description() ).c_str() );
+
+  checkRead(); requestAtoms(atom); 
+}
+
+void VolumeInSphere::setupRegions(){ }
+
+double VolumeInSphere::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
+  // Calculate position of atom wrt to origin
+  Vector fpos=pbcDistance( getPosition(0), cpos );
+  double dfunc, value = switchingFunction.calculateSqr( fpos.modulo2(), dfunc );    
+  derivatives.zero(); derivatives = dfunc*fpos; refders[0] = -derivatives;
+  // Add a virial contribution
+  vir -= Tensor(fpos,derivatives);
+  return value;
+}
+
+}
+}