Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2023 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 gridtools {
     110             : 
     111             : class FindSphericalContour : public ContourFindingBase {
     112             : private:
     113             :   unsigned nbins;
     114             :   double min, max;
     115             : public:
     116             :   static void registerKeywords( Keywords& keys );
     117             :   explicit FindSphericalContour(const ActionOptions&ao);
     118           4 :   unsigned getNumberOfQuantities() const override { return 2; }
     119             :   void compute( const unsigned& current, MultiValue& myvals ) const override;
     120             : };
     121             : 
     122       10421 : PLUMED_REGISTER_ACTION(FindSphericalContour,"FIND_SPHERICAL_CONTOUR")
     123             : 
     124           2 : void FindSphericalContour::registerKeywords( Keywords& keys ) {
     125           2 :   ContourFindingBase::registerKeywords( keys );
     126           4 :   keys.add("compulsory","NPOINTS","the number of points for which we are looking for the contour");
     127           4 :   keys.add("compulsory","INNER_RADIUS","the minimum radius on which to look for the contour");
     128           4 :   keys.add("compulsory","OUTER_RADIUS","the outer radius on which to look for the contour");
     129           4 :   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");
     130           2 : }
     131             : 
     132           1 : FindSphericalContour::FindSphericalContour(const ActionOptions&ao):
     133             :   Action(ao),
     134           1 :   ContourFindingBase(ao)
     135             : {
     136           1 :   if( ingrid->getDimension()!=3 ) error("input grid must be three dimensional");
     137             : 
     138           1 :   unsigned npoints; parse("NPOINTS",npoints);
     139           1 :   log.printf("  searching for %u points on dividing surface \n",npoints);
     140           3 :   parse("INNER_RADIUS",min); parse("OUTER_RADIUS",max); parse("NBINS",nbins);
     141           1 :   log.printf("  expecting to find dividing surface at radii between %f and %f \n",min,max);
     142           1 :   log.printf("  looking for contour in windows of length %f \n", (max-min)/nbins);
     143             :   // Set this here so the same set of grid points are used on every turn
     144           1 :   std::string vstring = "TYPE=fibonacci COMPONENTS=" + getLabel() + " COORDINATES=x,y,z PBC=F,F,F";
     145           2 :   auto grid=createGrid( "grid", vstring ); grid->setNoDerivatives();
     146           1 :   setAveragingAction( std::move(grid), true );
     147             :   // use mygrid, since at this point grid has been moved
     148           1 :   mygrid->setupFibonacciGrid( npoints );
     149             : 
     150           1 :   checkRead();
     151             :   // Create a task list
     152         101 :   for(unsigned i=0; i<npoints; ++i) addTaskToList( i );
     153           1 :   deactivateAllTasks();
     154         101 :   for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
     155           1 :   lockContributors();
     156           2 : }
     157             : 
     158         100 : void FindSphericalContour::compute( const unsigned& current, MultiValue& myvals ) const {
     159             :   // Generate contour point on inner sphere
     160         100 :   std::vector<double> contour_point(3), direction(3), der(3), tmp(3);
     161             :   // Retrieve this contour point from grid
     162         100 :   mygrid->getGridPointCoordinates( current, direction );
     163             :   // Now setup contour point on inner sphere
     164         400 :   for(unsigned j=0; j<3; ++j) {
     165         300 :     contour_point[j] = min*direction[j];
     166         300 :     direction[j] = (max-min)*direction[j] / static_cast<double>(nbins);
     167             :   }
     168             :   bool found=false;
     169         100 :   for(unsigned k=0; k<nbins; ++k) {
     170         400 :     for(unsigned j=0; j<3; ++j) tmp[j] = contour_point[j] + direction[j];
     171             :     double val1 = getDifferenceFromContour( contour_point, der );
     172             :     double val2 = getDifferenceFromContour( tmp, der );
     173         100 :     if( val1*val2<0 ) {
     174             :       findContour( direction, contour_point );
     175         400 :       double norm=0; for(unsigned j=0; j<3; ++j) norm += contour_point[j]*contour_point[j];
     176         100 :       myvals.setValue( 1, std::sqrt(norm) ); found=true; break;
     177             :     }
     178           0 :     for(unsigned j=0; j<3; ++j) contour_point[j] = tmp[j];
     179             :   }
     180           0 :   if( !found ) error("range does not bracket the dividing surface");
     181         100 : }
     182             : 
     183             : }
     184             : }