Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-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 "ContourFindingBase.h"
      24             : #include "tools/Random.h"
      25             : 
      26             : //+PLUMEDOC GRIDANALYSIS FIND_SPHERICAL_CONTOUR
      27             : /*
      28             : Find an isocontour in a three dimensional grid by searching over a Fibonacci sphere.
      29             : 
      30             : As discussed in the part of the manual on \ref Analysis PLUMED contains a number of tools that allow you to calculate
      31             : a function on a grid.  The function on this grid might be a \ref HISTOGRAM as a function of a few collective variables
      32             : or it might be a phase field that has been calculated using \ref MULTICOLVARDENS.  If this function has one or two input
      33             : arguments it is relatively straightforward to plot the function.  If by contrast the data has a three dimensions it can be
      34             : difficult to visualize.
      35             : 
      36             : This action provides one tool for visualizing these functions.  It can be used to search for a set of points on a contour
      37             : where the function takes a particular value.  In other words, for the function \f$f(x,y,z)\f$ this action would find a set
      38             : of points \f$\{x_c,y_c,z_c\}\f$ that have:
      39             : 
      40             : \f[
      41             : f(x_c,y_c,z_c) - c = 0
      42             : \f]
      43             : 
      44             : where \f$c\f$ is some constant value that is specified by the user.  The points on this contour are find by searching along a
      45             : set of equally spaced radii of a sphere that centered at on particular, user-specified atom or virtual atom.  To ensure that
      46             : these search radii are equally spaced on the surface of the sphere the search directions are generated by using a Fibonacci
      47             : spiral projected on a sphere.  In other words, the search directions are given by:
      48             : 
      49             : \f[
      50             : \mathbf{r}_i = \left(
      51             : \begin{matrix}
      52             : \sqrt{1 - y^2} \cos(\phi) \\
      53             : \frac{2i}{n} - 1 + \frac{1}{n}  \\
      54             : \sqrt{1 - y^2} \sin(\phi)
      55             : \end{matrix}
      56             : \right)
      57             : \f]
      58             : 
      59             : where \f$y\f$ is the quantity second component of the vector defined above, \f$n\f$ is the number of directions to look in and \f$\phi\f$ is
      60             : 
      61             : \f[
      62             : \phi = \mod(i + R, n) \pi ( 3 - \sqrt{5} )
      63             : \f]
      64             : 
      65             : where \f$R\f$ is a random variable between 0 and \f$n-1\f$ that is generated during the read in of the input file and that is fixed during
      66             : the whole calculation.
      67             : 
      68             : It is important to note that this action can only be used to detect contours in three dimensional functions.  In addition, this action will fail to
      69             : find the full set of contour  points if the contour does not have the same topology as a sphere.  If you are uncertain that the isocontours in your
      70             : function have a spherical topology you should use \ref FIND_CONTOUR in place of \ref FIND_SPHERICAL_CONTOUR.
      71             : 
      72             : \par Examples
      73             : 
      74             : The following input demonstrates how this action can be used.  The input here is used to study the shape of a droplet that has been formed during the
      75             : condensation of Lennard Jones from the vapor.  The input below achieves this by calculating the coordination numbers of all the atoms within the gas.
      76             : Obviously, those atoms within the droplet will have a large value for the coordination number while the isolated atoms in the gas will have a low value.
      77             : As such we can detect the sizes of the droplets by constructing a \ref CONTACT_MATRIX whose \f$ij\f$ element tells us whether atom \f$i\f$ and atom \f$j\f$
      78             : have coordination number that is greater that two.  The atoms within the various droplets within the system can then be found by performing a
      79             : \ref DFSCLUSTERING on this matrix to detect the connected components.  We can take the largest of these connected components and find the center of the droplet
      80             : by exploiting the functionality within \ref CENTER_OF_MULTICOLVAR.  We can then construct a phase field based on the positions of the atoms in the largest
      81             : cluster and the values of the coordination numbers of these atoms.  The final line in the input then finds the a set of points on the dividing surface that separates
      82             : the droplet from the surrounding gas.  The value of the phase field on this isocontour is equal to 0.75.
      83             : 
      84             : \plumedfile
      85             : # Calculate coordination numbers
      86             : c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
      87             : # Select coordination numbers that are more than 2.0
      88             : cf: MFILTER_MORE DATA=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
      89             : # Build a contact matrix
      90             : mat: CONTACT_MATRIX ATOMS=cf SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
      91             : # Find largest cluster
      92             : dfs: DFSCLUSTERING MATRIX=mat LOWMEM
      93             : clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1
      94             : # Find center of largest cluster
      95             : trans1: MTRANSFORM_MORE DATA=clust1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
      96             : cent: CENTER_OF_MULTICOLVAR DATA=trans1
      97             : # Calculate the phase field of the coordination
      98             : dens: MULTICOLVARDENS DATA=trans1 ORIGIN=cent DIR=xyz NBINS=30,30,30 BANDWIDTH=2.0,2.0,2.0
      99             : # Find the isocontour around the nucleus
     100             : sc: FIND_SPHERICAL_CONTOUR GRID=dens CONTOUR=0.85 INNER_RADIUS=10.0 OUTER_RADIUS=40.0 NPOINTS=100
     101             : # And print the grid to a file
     102             : GRID_TO_XYZ GRID=sc FILE=mysurface.xyz UNITS=A
     103             : \endplumedfile
     104             : 
     105             : */
     106             : //+ENDPLUMEDOC
     107             : 
     108             : namespace PLMD {
     109             : namespace contour {
     110             : 
     111             : class FindSphericalContour : public ContourFindingBase {
     112             : private:
     113             :   unsigned nbins, npoints;
     114             :   double min, max;
     115             :   gridtools::GridCoordinatesObject gridcoords;
     116             : public:
     117             :   static void registerKeywords( Keywords& keys );
     118             :   explicit FindSphericalContour(const ActionOptions&ao);
     119           1 :   void setupValuesOnFirstStep() override {}
     120             :   unsigned getNumberOfDerivatives() override ;
     121             :   std::vector<std::string> getGridCoordinateNames() const override ;
     122             :   const gridtools::GridCoordinatesObject& getGridCoordinatesObject() const override ;
     123             :   void performTask( const unsigned& current, MultiValue& myvals ) const override;
     124             :   void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     125             :                           const unsigned& bufstart, std::vector<double>& buffer ) const override ;
     126             : };
     127             : 
     128             : PLUMED_REGISTER_ACTION(FindSphericalContour,"FIND_SPHERICAL_CONTOUR")
     129             : 
     130           3 : void FindSphericalContour::registerKeywords( Keywords& keys ) {
     131           3 :   ContourFindingBase::registerKeywords( keys );
     132           3 :   keys.add("compulsory","NPOINTS","the number of points for which we are looking for the contour");
     133           3 :   keys.add("compulsory","INNER_RADIUS","the minimum radius on which to look for the contour");
     134           3 :   keys.add("compulsory","OUTER_RADIUS","the outer radius on which to look for the contour");
     135           3 :   keys.add("compulsory","NBINS","1","the number of discrete sections in which to divide the distance between the inner and outer radius when searching for a contour");
     136           6 :   keys.setValueDescription("grid","a grid on a Fibonacci sphere that describes the radial distance from the origin for the points on the Willard-Chandler surface");
     137           3 : }
     138             : 
     139           1 : FindSphericalContour::FindSphericalContour(const ActionOptions&ao):
     140             :   Action(ao),
     141           1 :   ContourFindingBase(ao) {
     142           1 :   if( getPntrToArgument(0)->getRank()!=3 ) {
     143           0 :     error("input grid must be three dimensional");
     144             :   }
     145             : 
     146           1 :   parse("NPOINTS",npoints);
     147           1 :   log.printf("  searching for %u points on dividing surface \n",npoints);
     148           1 :   parse("INNER_RADIUS",min);
     149           1 :   parse("OUTER_RADIUS",max);
     150           1 :   parse("NBINS",nbins);
     151           1 :   log.printf("  expecting to find dividing surface at radii between %f and %f \n",min,max);
     152           1 :   log.printf("  looking for contour in windows of length %f \n", (max-min)/nbins);
     153             :   // Set this here so the same set of grid points are used on every turn
     154           1 :   std::vector<bool> ipbc( 3, false );
     155           1 :   gridcoords.setup( "fibonacci", ipbc, npoints, 0.0 );
     156             :   // Now create a value
     157           1 :   std::vector<unsigned> shape( 3 );
     158           1 :   shape[0]=npoints;
     159           1 :   shape[1]=shape[2]=1;
     160           1 :   addValueWithDerivatives( shape );
     161           1 :   setNotPeriodic();
     162           1 :   getPntrToComponent(0)->buildDataStore();
     163           1 :   checkRead();
     164           1 : }
     165             : 
     166           2 : unsigned FindSphericalContour::getNumberOfDerivatives() {
     167           2 :   return gridcoords.getDimension();
     168             : }
     169             : 
     170           0 : std::vector<std::string> FindSphericalContour::getGridCoordinateNames() const {
     171           0 :   gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
     172           0 :   plumed_assert( ag );
     173           0 :   return ag->getGridCoordinateNames();
     174             : }
     175             : 
     176           3 : const gridtools::GridCoordinatesObject& FindSphericalContour::getGridCoordinatesObject() const {
     177           3 :   return gridcoords;
     178             : }
     179             : 
     180         200 : void FindSphericalContour::performTask( const unsigned& current, MultiValue& myvals ) const {
     181             :   // Generate contour point on inner sphere
     182         200 :   std::vector<double> contour_point(3), direction(3), der(3), tmp(3);
     183             :   // Retrieve this contour point from grid
     184         200 :   gridcoords.getGridPointCoordinates( current, direction );
     185             :   // Now setup contour point on inner sphere
     186         800 :   for(unsigned j=0; j<3; ++j) {
     187         600 :     contour_point[j] = min*direction[j];
     188         600 :     direction[j] = (max-min)*direction[j] / static_cast<double>(nbins);
     189             :   }
     190             : 
     191             :   bool found=false;
     192         200 :   for(unsigned k=0; k<nbins; ++k) {
     193         800 :     for(unsigned j=0; j<3; ++j) {
     194         600 :       tmp[j] = contour_point[j] + direction[j];
     195             :     }
     196         200 :     double val1 = getDifferenceFromContour( contour_point, der );
     197         200 :     double val2 = getDifferenceFromContour( tmp, der );
     198         200 :     if( val1*val2<0 ) {
     199             :       findContour( direction, contour_point );
     200             :       double norm=0;
     201         800 :       for(unsigned j=0; j<3; ++j) {
     202         600 :         norm += contour_point[j]*contour_point[j];
     203             :       }
     204         200 :       myvals.setValue( getConstPntrToComponent(0)->getPositionInStream(), sqrt(norm) );
     205             :       found=true;
     206             :       break;
     207             :     }
     208           0 :     for(unsigned j=0; j<3; ++j) {
     209           0 :       contour_point[j] = tmp[j];
     210             :     }
     211             :   }
     212             :   if( !found ) {
     213           0 :     error("range does not bracket the dividing surface");
     214             :   }
     215         200 : }
     216             : 
     217         200 : void FindSphericalContour::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     218             :     const unsigned& bufstart, std::vector<double>& buffer ) const {
     219         200 :   plumed_assert( valindex==0 );
     220         200 :   unsigned istart = bufstart + (1+gridcoords.getDimension())*code;
     221         200 :   unsigned valout = getConstPntrToComponent(0)->getPositionInStream();
     222         200 :   buffer[istart] += myvals.get( valout );
     223         200 : }
     224             : 
     225             : }
     226             : }