Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "ActionVolume.h"
      26             : #include "VolumeShortcut.h"
      27             : 
      28             : //+PLUMEDOC VOLUMES INSPHERE
      29             : /*
      30             : 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.
      31             : 
      32             : This action can be used to calculate whether each of the atoms are within a particular part of the simulation box or not as illustrated by the following example:
      33             : 
      34             : ```plumed
      35             : f: FIXEDATOM AT=0,0,0
      36             : a: INSPHERE ATOMS=1-100 CENTER=f RADIUS={GAUSSIAN R_0=1.5}
      37             : PRINT ARG=a FILE=colvar
      38             : ```
      39             : 
      40             : The 100 elements of the vector `a` that is returned from the INSPHERE action in the above input are calculated using:
      41             : 
      42             : $$
      43             : w(x_i,y_i,z_i) = \sigma\left( \sqrt{ x_i^2 + y_i^2 + z_i^2}  \right)
      44             : $$
      45             : 
      46             : In this expression $x_i$, $y_i$ and $z_i$ are the components of the vector that connects the $i$th atom that was specified using the `ATOMS` keyword to the atom that was specified using the `CENTER` keyword and
      47             : $\sigma$ is one of the switching functions that is described that in the documentation for the action [LESS_THAN](LESS_THAN.md).  In other words,
      48             : $w(x_i,y_i,z_i)$ is 1 if atom $i$ is within a sphere with a radial extent that is determined by the parameters of the specified switching function
      49             : and zero otherwise.
      50             : 
      51             : If you want to caluclate and print the number of atoms in the sphere of interest you can use an input like the one shown below:
      52             : 
      53             : ```plumed
      54             : f: FIXEDATOM AT=0,0,0
      55             : a: INSPHERE ATOMS=1-100 CENTER=f RADIUS={GAUSSIAN R_0=1.5}
      56             : s: SUM ARG=a PERIODIC=NO
      57             : PRINT ARG=s FILE=colvar
      58             : ```
      59             : 
      60             : You can also calculate the average values of symmetry functions in the sphere of interest by using inputs similar to those described the documentation for the [AROUND](AROUND.md)
      61             : action. In other words, you can swap out AROUND actions for an INSPHERE actions.  Also as with [AROUND](AROUND.md), earlier versions of PLUMED used a different syntax for doing these types of calculations, which can
      62             : still be used with this new version of the command.  However, we strongly recommend using the newer syntax.
      63             : 
      64             : 
      65             : */
      66             : //+ENDPLUMEDOC
      67             : 
      68             : //+PLUMEDOC MCOLVAR INSPHERE_CALC
      69             : /*
      70             : Calculate a vector from the input positions with elements equal to one when the positions are in a particular part of the cell and elements equal to zero otherwise
      71             : 
      72             : \par Examples
      73             : 
      74             : */
      75             : //+ENDPLUMEDOC
      76             : 
      77             : namespace PLMD {
      78             : namespace volumes {
      79             : 
      80             : class VolumeInSphere : public ActionVolume {
      81             : private:
      82             :   Vector origin;
      83             :   SwitchingFunction switchingFunction;
      84             : public:
      85             :   static void registerKeywords( Keywords& keys );
      86             :   explicit VolumeInSphere(const ActionOptions& ao);
      87             :   void setupRegions() override;
      88             :   double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
      89             : };
      90             : 
      91             : PLUMED_REGISTER_ACTION(VolumeInSphere,"INSPHERE_CALC")
      92             : char glob_sphere[] = "INSPHERE";
      93             : typedef VolumeShortcut<glob_sphere> VolumeInSphereShortcut;
      94             : PLUMED_REGISTER_ACTION(VolumeInSphereShortcut,"INSPHERE")
      95             : 
      96          36 : void VolumeInSphere::registerKeywords( Keywords& keys ) {
      97          36 :   ActionVolume::registerKeywords( keys );
      98          72 :   keys.setDisplayName("INSPHERE");
      99          36 :   keys.add("atoms","CENTER","the atom whose vicinity we are interested in examining");
     100          36 :   keys.add("atoms-2","ATOM","the atom whose vicinity we are interested in examining");
     101          36 :   keys.add("compulsory","RADIUS","the switching function that tells us the extent of the sphereical region of interest");
     102          36 :   keys.remove("SIGMA");
     103          72 :   keys.linkActionInDocs("RADIUS","LESS_THAN");
     104          36 : }
     105             : 
     106          10 : VolumeInSphere::VolumeInSphere(const ActionOptions& ao):
     107             :   Action(ao),
     108          10 :   ActionVolume(ao) {
     109             :   std::vector<AtomNumber> atom;
     110          20 :   parseAtomList("CENTER",atom);
     111          10 :   if( atom.size()==0 ) {
     112           0 :     parseAtomList("ATOM",atom);
     113             :   }
     114          10 :   if( atom.size()!=1 ) {
     115           0 :     error("should only be one atom specified");
     116             :   }
     117          10 :   log.printf("  center of sphere is at position of atom : %d\n",atom[0].serial() );
     118             : 
     119             :   std::string sw, errors;
     120          20 :   parse("RADIUS",sw);
     121          10 :   if(sw.length()==0) {
     122           0 :     error("missing RADIUS keyword");
     123             :   }
     124          10 :   switchingFunction.set(sw,errors);
     125          10 :   if( errors.length()!=0 ) {
     126           0 :     error("problem reading RADIUS keyword : " + errors );
     127             :   }
     128          10 :   log.printf("  radius of sphere is given by %s \n", ( switchingFunction.description() ).c_str() );
     129             : 
     130          10 :   checkRead();
     131          10 :   requestAtoms(atom);
     132          10 : }
     133             : 
     134         120 : void VolumeInSphere::setupRegions() { }
     135             : 
     136      423194 : double VolumeInSphere::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
     137             :   // Calculate position of atom wrt to origin
     138      423194 :   Vector fpos=pbcDistance( getPosition(0), cpos );
     139      423194 :   double dfunc, value = switchingFunction.calculateSqr( fpos.modulo2(), dfunc );
     140      423194 :   derivatives.zero();
     141      423194 :   derivatives = dfunc*fpos;
     142      423194 :   refders[0] = -derivatives;
     143             :   // Add a virial contribution
     144      423194 :   vir -= Tensor(fpos,derivatives);
     145      423194 :   return value;
     146             : }
     147             : 
     148             : }
     149             : }