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           6 :   keys.add("compulsory","NPOINTS","the number of points for which we are looking for the contour");
     133           6 :   keys.add("compulsory","INNER_RADIUS","the minimum radius on which to look for the contour");
     134           6 :   keys.add("compulsory","OUTER_RADIUS","the outer radius on which to look for the contour");
     135           6 :   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           3 :   keys.setValueDescription("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             : {
     143           1 :   if( getPntrToArgument(0)->getRank()!=3 ) error("input grid must be three dimensional");
     144             : 
     145           2 :   parse("NPOINTS",npoints); log.printf("  searching for %u points on dividing surface \n",npoints);
     146           3 :   parse("INNER_RADIUS",min); parse("OUTER_RADIUS",max); parse("NBINS",nbins);
     147           1 :   log.printf("  expecting to find dividing surface at radii between %f and %f \n",min,max);
     148           1 :   log.printf("  looking for contour in windows of length %f \n", (max-min)/nbins);
     149             :   // Set this here so the same set of grid points are used on every turn
     150           1 :   std::vector<bool> ipbc( 3, false ); gridcoords.setup( "fibonacci", ipbc, npoints, 0.0 );
     151             :   // Now create a value
     152           1 :   std::vector<unsigned> shape( 3 ); shape[0]=npoints; shape[1]=shape[2]=1;
     153           1 :   addValueWithDerivatives( shape ); setNotPeriodic();
     154           1 :   getPntrToComponent(0)->buildDataStore(); checkRead();
     155           1 : }
     156             : 
     157           2 : unsigned FindSphericalContour::getNumberOfDerivatives() {
     158           2 :   return gridcoords.getDimension();
     159             : }
     160             : 
     161           0 : std::vector<std::string> FindSphericalContour::getGridCoordinateNames() const {
     162           0 :   gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
     163           0 :   plumed_assert( ag ); return ag->getGridCoordinateNames();
     164             : }
     165             : 
     166           3 : const gridtools::GridCoordinatesObject& FindSphericalContour::getGridCoordinatesObject() const {
     167           3 :   return gridcoords;
     168             : }
     169             : 
     170         200 : void FindSphericalContour::performTask( const unsigned& current, MultiValue& myvals ) const {
     171             :   // Generate contour point on inner sphere
     172         200 :   std::vector<double> contour_point(3), direction(3), der(3), tmp(3);
     173             :   // Retrieve this contour point from grid
     174         200 :   gridcoords.getGridPointCoordinates( current, direction );
     175             :   // Now setup contour point on inner sphere
     176         800 :   for(unsigned j=0; j<3; ++j) {
     177         600 :     contour_point[j] = min*direction[j];
     178         600 :     direction[j] = (max-min)*direction[j] / static_cast<double>(nbins);
     179             :   }
     180             : 
     181             :   bool found=false;
     182         200 :   for(unsigned k=0; k<nbins; ++k) {
     183         800 :     for(unsigned j=0; j<3; ++j) tmp[j] = contour_point[j] + direction[j];
     184         200 :     double val1 = getDifferenceFromContour( contour_point, der );
     185         200 :     double val2 = getDifferenceFromContour( tmp, der );
     186         200 :     if( val1*val2<0 ) {
     187             :       findContour( direction, contour_point );
     188         800 :       double norm=0; for(unsigned j=0; j<3; ++j) norm += contour_point[j]*contour_point[j];
     189         200 :       myvals.setValue( getConstPntrToComponent(0)->getPositionInStream(), sqrt(norm) ); found=true; break;
     190             :     }
     191           0 :     for(unsigned j=0; j<3; ++j) contour_point[j] = tmp[j];
     192             :   }
     193           0 :   if( !found ) error("range does not bracket the dividing surface");
     194         200 : }
     195             : 
     196         200 : void FindSphericalContour::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     197             :     const unsigned& bufstart, std::vector<double>& buffer ) const {
     198         200 :   plumed_assert( valindex==0 ); unsigned istart = bufstart + (1+gridcoords.getDimension())*code;
     199         200 :   unsigned valout = getConstPntrToComponent(0)->getPositionInStream(); buffer[istart] += myvals.get( valout );
     200         200 : }
     201             : 
     202             : }
     203             : }