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             : #ifdef __PLUMED_HAS_OPENACC
      23             : #define __PLUMED_USE_OPENACC 1
      24             : #endif //__PLUMED_HAS_OPENACC
      25             : #include "core/ActionRegister.h"
      26             : #include "tools/Pbc.h"
      27             : #include "tools/SwitchingFunction.h"
      28             : #include "ActionVolume.h"
      29             : #include "VolumeShortcut.h"
      30             : 
      31             : //+PLUMEDOC VOLUMES INSPHERE
      32             : /*
      33             : 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.
      34             : 
      35             : 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:
      36             : 
      37             : ```plumed
      38             : f: FIXEDATOM AT=0,0,0
      39             : a: INSPHERE ATOMS=1-100 CENTER=f RADIUS={GAUSSIAN R_0=1.5}
      40             : PRINT ARG=a FILE=colvar
      41             : ```
      42             : 
      43             : The 100 elements of the vector `a` that is returned from the INSPHERE action in the above input are calculated using:
      44             : 
      45             : $$
      46             : w(x_i,y_i,z_i) = \sigma\left( \sqrt{ x_i^2 + y_i^2 + z_i^2}  \right)
      47             : $$
      48             : 
      49             : 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
      50             : $\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,
      51             : $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
      52             : and zero otherwise.
      53             : 
      54             : 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:
      55             : 
      56             : ```plumed
      57             : f: FIXEDATOM AT=0,0,0
      58             : a: INSPHERE ATOMS=1-100 CENTER=f RADIUS={GAUSSIAN R_0=1.5}
      59             : s: SUM ARG=a PERIODIC=NO
      60             : PRINT ARG=s FILE=colvar
      61             : ```
      62             : 
      63             : If by constrast you want to calculate and print the number of atoms that are not in the sphere of interest you OUTSIDE flag as shown below
      64             : 
      65             : ```plumed
      66             : f: FIXEDATOM AT=0,0,0
      67             : a: INSPHERE ...
      68             :   ATOMS=1-100 CENTER=f
      69             :   RADIUS={GAUSSIAN R_0=1.5}
      70             :   OUTSIDE
      71             : ...
      72             : s: SUM ARG=a PERIODIC=NO
      73             : PRINT ARG=s FILE=colvar
      74             : ```
      75             : 
      76             : !!! note ""
      77             : 
      78             :     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)
      79             :     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
      80             :     still be used with this new version of the command.  We strongly recommend using the newer syntax but if you are interested in the
      81             :     old syntax you can find more information in the old syntax section of the documentation for [AROUND](AROUND.md).
      82             : 
      83             : */
      84             : //+ENDPLUMEDOC
      85             : 
      86             : namespace PLMD {
      87             : namespace volumes {
      88             : 
      89          60 : struct VolumeInSphere {
      90             : #ifdef __PLUMED_USE_OPENACC
      91             :   SwitchingFunctionAccelerable switchingFunction;
      92             : #else
      93             :   SwitchingFunction switchingFunction;
      94             : #endif //__PLUMED_USE_OPENACC
      95             :   static void registerKeywords( Keywords& keys );
      96             :   void parseInput( ActionVolume<VolumeInSphere>* action );
      97             :   void setupRegions( ActionVolume<VolumeInSphere>* action, const Pbc& pbc, const std::vector<Vector>& positions ) {}
      98             :   static void parseAtoms( ActionVolume<VolumeInSphere>* action, std::vector<AtomNumber>& atom );
      99             :   static void calculateNumberInside( const VolumeInput& input, const VolumeInSphere& actioninput, VolumeOutput& output );
     100             : #ifdef __PLUMED_USE_OPENACC
     101             :   void toACCDevice() const {
     102             : #pragma acc enter data copyin(this[0:1])
     103             :     switchingFunction.toACCDevice();
     104             :   }
     105             :   void removeFromACCDevice() const {
     106             :     switchingFunction.removeFromACCDevice();
     107             : #pragma acc exit data delete(this[0:1])
     108             :   }
     109             : #endif //__PLUMED_USE_OPENACC
     110             : };
     111             : 
     112             : typedef ActionVolume<VolumeInSphere> Vols;
     113             : PLUMED_REGISTER_ACTION(Vols,"INSPHERE_CALC")
     114             : char glob_sphere[] = "INSPHERE";
     115             : typedef VolumeShortcut<glob_sphere> VolumeInSphereShortcut;
     116             : PLUMED_REGISTER_ACTION(VolumeInSphereShortcut,"INSPHERE")
     117             : 
     118          51 : void VolumeInSphere::registerKeywords( Keywords& keys ) {
     119         102 :   keys.setDisplayName("INSPHERE");
     120          51 :   keys.add("atoms","CENTER","the atom whose vicinity we are interested in examining");
     121          51 :   keys.addDeprecatedKeyword("ATOM","CENTER");
     122          51 :   keys.add("compulsory","RADIUS","the switching function that tells us the extent of the sphereical region of interest");
     123         102 :   keys.linkActionInDocs("RADIUS","LESS_THAN");
     124          51 : }
     125             : 
     126          15 : void VolumeInSphere::parseInput( ActionVolume<VolumeInSphere>* action ) {
     127             :   std::string errors;
     128             :   std::string swinput;
     129          30 :   action->parse("RADIUS",swinput);
     130          15 :   if(swinput.length()==0) {
     131           0 :     action->error("missing RADIUS keyword");
     132             :   }
     133             : 
     134          15 :   switchingFunction.set(swinput,errors);
     135          15 :   if( errors.length()!=0 ) {
     136           0 :     action->error("problem reading RADIUS keyword : " + errors );
     137             :   }
     138             : 
     139          15 :   action->log.printf("  radius of sphere is given by %s \n",
     140          30 :                      switchingFunction.description().c_str() );
     141          15 : }
     142             : 
     143          15 : void VolumeInSphere::parseAtoms( ActionVolume<VolumeInSphere>* action, std::vector<AtomNumber>& atom ) {
     144          30 :   action->parseAtomList("CENTER",atom);
     145          15 :   if( atom.size()==0 ) {
     146           0 :     action->parseAtomList("ATOM",atom);
     147             :   }
     148          15 :   if( atom.size()!=1 ) {
     149           0 :     action->error("should only be one atom specified");
     150             :   }
     151          15 :   action->log.printf("  center of sphere is at position of atom : %d\n",atom[0].serial() );
     152          15 : }
     153             : 
     154      946046 : void VolumeInSphere::calculateNumberInside( const VolumeInput& input,
     155             :     const VolumeInSphere& actioninput,
     156             :     VolumeOutput& output ) {
     157             :   // Calculate position of atom wrt to origin
     158      946046 :   Vector fpos=input.pbc.distance( Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]), Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
     159             :   double dfunc;
     160      946046 :   output.values[0] = actioninput.switchingFunction.calculateSqr( fpos.modulo2(), dfunc );
     161      946046 :   output.derivatives = dfunc*fpos;
     162      946046 :   output.refders[0][0] = -output.derivatives[0];
     163      946046 :   output.refders[0][1] = -output.derivatives[1];
     164      946046 :   output.refders[0][2] = -output.derivatives[2];
     165             :   // Add a virial contribution
     166     1892092 :   output.virial.set( 0, -Tensor(fpos,Vector(output.derivatives[0], output.derivatives[1], output.derivatives[2])) );
     167      946046 : }
     168             : 
     169             : }
     170             : }