Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2017-2020 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 "core/ActionRegister.h"
      23             : #include "tools/Pbc.h"
      24             : #include "tools/SwitchingFunction.h"
      25             : #include "tools/LinkCells.h"
      26             : #include "ActionVolume.h"
      27             : #include "VolumeShortcut.h"
      28             : #include <memory>
      29             : 
      30             : //+PLUMEDOC VOLUMES INENVELOPE
      31             : /*
      32             : This quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a region where the density of a certain type of atom is high.
      33             : 
      34             : This collective variable can be used to determine whether colvars are within region where the density
      35             : of a particular atom is high.  This is achieved by calculating the following function at the point where
      36             : the atom is located \f$(x,y,z)\f$:
      37             : 
      38             : \f[
      39             : w_j = 1 - \sigma\left[ \sum_{i=1}^N K\left( \frac{x-x_i}{\sigma_x},\frac{y-y_i}{\sigma_y},\frac{z-z_i}{\sigma_z} \right) \right]
      40             : \f]
      41             : 
      42             : Here \f$\sigma\f$ is a \ref switchingfunction and \f$K\f$ is a \ref kernelfunctions.  The sum runs over the atoms
      43             : specified using the ATOMS keyword and a \f$w_j\f$ value is calculated for each of the central atoms of the input
      44             : multicolvar.
      45             : 
      46             : \par Examples
      47             : 
      48             : The input below calculates a density field from the positions of atoms 1-14400.  The number of the atoms
      49             : that are specified in the DENSITY action that are within a region where the density field is greater than
      50             : 2.0 is then calculated.
      51             : 
      52             : \plumedfile
      53             : d1: DENSITY SPECIES=14401-74134:3 LOWMEM
      54             : fi: INENVELOPE DATA=d1 ATOMS=1-14400 CONTOUR={RATIONAL D_0=2.0 R_0=1.0} BANDWIDTH=0.1,0.1,0.1 LOWMEM
      55             : PRINT ARG=fi FILE=colvar
      56             : \endplumedfile
      57             : 
      58             : */
      59             : //+ENDPLUMEDOC
      60             : 
      61             : //+PLUMEDOC VOLUMES INENVELOPE_CALC
      62             : /*
      63             : This quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a region where the density of a certain type of atom is high.
      64             : 
      65             : \par Examples
      66             : 
      67             : */
      68             : //+ENDPLUMEDOC
      69             : 
      70             : namespace PLMD {
      71             : namespace volumes {
      72             : 
      73             : class VolumeInEnvelope : public ActionVolume {
      74             : private:
      75             :   LinkCells mylinks;
      76             :   double gvol;
      77             :   std::vector<std::unique_ptr<Value>> pos;
      78             :   std::vector<Vector> ltmp_pos;
      79             :   std::vector<unsigned> ltmp_ind;
      80             :   std::vector<double> bandwidth;
      81             :   SwitchingFunction sfunc, switchingFunction;
      82             : public:
      83             :   static void registerKeywords( Keywords& keys );
      84             :   explicit VolumeInEnvelope(const ActionOptions& ao);
      85             :   void setupRegions() override;
      86             :   double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
      87             : };
      88             : 
      89             : PLUMED_REGISTER_ACTION(VolumeInEnvelope,"INENVELOPE_CALC")
      90             : char glob_contours[] = "INENVELOPE";
      91             : typedef VolumeShortcut<glob_contours> VolumeInEnvelopeShortcut;
      92             : PLUMED_REGISTER_ACTION(VolumeInEnvelopeShortcut,"INENVELOPE")
      93             : 
      94           7 : void VolumeInEnvelope::registerKeywords( Keywords& keys ) {
      95          14 :   ActionVolume::registerKeywords( keys ); keys.remove("SIGMA"); keys.setDisplayName("INENVELOPE");
      96          14 :   keys.add("atoms","FIELD_ATOMS","the atom whose positions we are constructing a field from");
      97          14 :   keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
      98          14 :   keys.add("compulsory","CONTOUR","a switching funciton that tells PLUMED how large the density should be");
      99          14 :   keys.add("compulsory","CUTOFF","6.25","the cutoff at which to stop evaluating the kernel functions is set equal to sqrt(2*x)*bandwidth in each direction where x is this number");
     100           7 : }
     101             : 
     102           1 : VolumeInEnvelope::VolumeInEnvelope(const ActionOptions& ao):
     103             :   Action(ao),
     104             :   ActionVolume(ao),
     105           1 :   mylinks(comm)
     106             : {
     107           1 :   std::vector<AtomNumber> atoms; parseAtomList("FIELD_ATOMS",atoms);
     108           1 :   log.printf("  creating density field from atoms : ");
     109           9 :   for(unsigned i=0; i<atoms.size(); ++i) log.printf("%d ",atoms[i].serial() );
     110           1 :   log.printf("\n"); ltmp_ind.resize( atoms.size() ); ltmp_pos.resize( atoms.size() );
     111           9 :   for(unsigned i=0; i<atoms.size(); ++i) ltmp_ind[i]=i;
     112             : 
     113           2 :   std::string sw, errors; parse("CONTOUR",sw);
     114           1 :   if(sw.length()==0) error("missing CONTOUR keyword");
     115           1 :   sfunc.set(sw,errors);
     116           1 :   if( errors.length()!=0 ) error("problem reading RADIUS keyword : " + errors );
     117           1 :   log.printf("  density at atom must be larger than %s \n", ( sfunc.description() ).c_str() );
     118             : 
     119           1 :   std::vector<double> pp(3,0.0); bandwidth.resize(3); parseVector("BANDWIDTH",bandwidth);
     120           2 :   log.printf("  using %s kernel with bandwidths %f %f %f \n",getKernelType().c_str(),bandwidth[0],bandwidth[1],bandwidth[2] );
     121           2 :   std::string errors2; switchingFunction.set("GAUSSIAN R_0=1.0 NOSTRETCH", errors2 );
     122           1 :   if( errors2.length()!=0 ) error("problem reading switching function description " + errors2);
     123           4 :   double det=1; for(unsigned i=0; i<bandwidth.size(); ++i) det*=bandwidth[i]*bandwidth[i];
     124           2 :   gvol=1.0; if( getKernelType()=="gaussian" ) gvol=pow( 2*pi, 0.5*bandwidth.size() ) * pow( det, 0.5 );
     125           0 :   else error("cannot use kernel other than gaussian");
     126           2 :   double dp2cutoff; parse("CUTOFF",dp2cutoff); double maxs =  sqrt(2*dp2cutoff)*bandwidth[0];
     127           3 :   for(unsigned j=1; j<bandwidth.size(); ++j) {
     128           2 :     if( sqrt(2*dp2cutoff)*bandwidth[j]>maxs ) maxs=sqrt(2*dp2cutoff)*bandwidth[j];
     129             :   }
     130           1 :   checkRead(); requestAtoms(atoms); mylinks.setCutoff( maxs );
     131           1 : }
     132             : 
     133           5 : void VolumeInEnvelope::setupRegions() {
     134          45 :   for(unsigned i=0; i<ltmp_ind.size(); ++i) { ltmp_pos[i]=getPosition(i); }
     135           5 :   mylinks.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
     136           5 : }
     137             : 
     138        1000 : double VolumeInEnvelope::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
     139        1000 :   unsigned ncells_required=0, natoms=1; std::vector<unsigned> cells_required( mylinks.getNumberOfCells() ), indices( 1 + getNumberOfAtoms() );
     140        1000 :   mylinks.addRequiredCells( mylinks.findMyCell( cpos ), ncells_required, cells_required );
     141        1000 :   indices[0]=getNumberOfAtoms(); mylinks.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
     142        1000 :   double value=0; std::vector<double> der(3); Vector tder;
     143             : 
     144             :   // convert pointer once
     145        1000 :   auto pos_ptr=Tools::unique2raw(pos);
     146        9000 :   for(unsigned i=1; i<natoms; ++i) {
     147        8000 :     Vector dist = pbcDistance( cpos, getPosition( indices[i] ) );
     148       32000 :     double dval=0; for(unsigned j=0; j<3; ++j) { der[j] = dist[j]/bandwidth[j]; dval += der[j]*der[j]; der[j] = der[j] / bandwidth[j]; }
     149        8000 :     double dfunc; value += switchingFunction.calculateSqr( dval, dfunc ) / gvol; double tmp = dfunc / gvol;
     150       32000 :     for(unsigned j=0; j<3; ++j) {
     151       24000 :       derivatives[j] -= tmp*der[j]; refders[ indices[i] ][j] += tmp*der[j]; tder[j]=tmp*der[j];
     152             :     }
     153        8000 :     vir -= Tensor( tder, dist );
     154             :   }
     155        1000 :   double deriv, fval = sfunc.calculate( value, deriv );
     156        1000 :   derivatives *= -deriv*value; vir *= -deriv*value;
     157        9000 :   for(unsigned i=1; i<natoms; ++i) refders[ indices[i] ] *= -deriv*value;
     158        2000 :   return 1.0 - fval;
     159             : }
     160             : 
     161             : }
     162             : }