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 "core/PlumedMain.h"
      25             : 
      26             : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR_SURFACE
      27             : /*
      28             : Find an isocontour by searching along either the x, y or z direction.
      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 lines
      45             : that run parallel to the \f$x\f$, \f$y\f$ or \f$z\f$ axis of the simulation cell.  The result is, therefore, a two dimensional
      46             : function evaluated on a grid that gives us the height of the interface as a function of two coordinates.
      47             : 
      48             : 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
      49             : find the full set of contour  points if the contour does not have the same topology as an infinite plane.  If you are uncertain that the isocontours in your
      50             : function have the appropriate topology you should use \ref FIND_CONTOUR in place of \ref FIND_CONTOUR_SURFACE.
      51             : 
      52             : 
      53             : \par Examples
      54             : 
      55             : The input shown below was used to analyze the results from a simulation of an interface between solid and molten Lennard Jones.  The interface between
      56             : the solid and the liquid was set up in the plane perpendicular to the \f$z\f$ direction of the simulation cell.   The input below calculates something
      57             : akin to a Willard-Chandler dividing surface \cite wcsurface between the solid phase and the liquid phase.  There are two of these interfaces within the
      58             : simulation box because of the periodic boundary conditions but we were able to determine that one of these two surfaces lies in a particular part of the
      59             : simulation box.  The input below detects the height profile of one of these two interfaces.  It does so by computing a phase field average of the
      60             : \ref FCCUBIC symmetry function using the \ref MULTICOLVARDENS action.  Notice that we use the fact that we know roughly where the interface is when specifying
      61             : how this phase field is to be calculated and specify the region over the \f$z\f$-axis in which we are going to search for the phase field in the line defining
      62             : the \ref MULTICOLVARDENS.  Once we have calculated the phase field we search for contour points on the lines that run parallel to the \f$z\f$-direction of the cell
      63             : box using the FIND_CONTOUR_SURFACE command.  The final result is a \f$14 \times 14\f$ grid of values for the height of the interface as a function of the \f$(x,y)\f$
      64             : position.  This grid is then output to a file called contour2.dat.
      65             : 
      66             : Notice that the commands below calculate the instantaneous position of the surface separating the solid and liquid and that as such the accumulated average is cleared
      67             : on every step.
      68             : 
      69             : \plumedfile
      70             : UNITS NATURAL
      71             : FCCUBIC ...
      72             :   SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
      73             :   ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619 LABEL=fcc
      74             : ... FCCUBIC
      75             : 
      76             : dens2: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,50 ZREDUCED ZLOWER=6.0 ZUPPER=11.0 BANDWIDTH=1.0,1.0,1.0 CLEAR=1
      77             : 
      78             : ss2: FIND_CONTOUR_SURFACE GRID=dens2 CONTOUR=0.42 SEARCHDIR=z STRIDE=1 CLEAR=1
      79             : DUMPGRID GRID=ss2 FILE=contour2.dat FMT=%8.4f STRIDE=1
      80             : \endplumedfile
      81             : 
      82             : */
      83             : //+ENDPLUMEDOC
      84             : 
      85             : namespace PLMD {
      86             : namespace contour {
      87             : 
      88             : class FindContourSurface : public ContourFindingBase {
      89             : private:
      90             :   unsigned dir_n;
      91             :   std::vector<unsigned> ones;
      92             :   std::vector<unsigned> gdirs;
      93             :   std::vector<double> direction;
      94             :   std::vector<std::string> gnames;
      95             :   gridtools::GridCoordinatesObject gridcoords;
      96             : public:
      97             :   static void registerKeywords( Keywords& keys );
      98             :   explicit FindContourSurface(const ActionOptions&ao);
      99             :   void setupValuesOnFirstStep() override;
     100             :   unsigned getNumberOfDerivatives() override ;
     101             :   std::vector<std::string> getGridCoordinateNames() const override ;
     102             :   const gridtools::GridCoordinatesObject& getGridCoordinatesObject() const override ;
     103             :   void performTask( const unsigned& current, MultiValue& myvals ) const override;
     104             :   void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     105             :                           const unsigned& bufstart, std::vector<double>& buffer ) const override ;
     106             : };
     107             : 
     108             : PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
     109             : 
     110           3 : void FindContourSurface::registerKeywords( Keywords& keys ) {
     111           3 :   ContourFindingBase::registerKeywords( keys );
     112           3 :   keys.add("compulsory","SEARCHDIR","In which directions do you wish to search for the contour.");
     113           6 :   keys.setValueDescription("grid","a grid containing the location of the points in the Willard-Chandler surface along the chosen direction");
     114           3 : }
     115             : 
     116           1 : FindContourSurface::FindContourSurface(const ActionOptions&ao):
     117             :   Action(ao),
     118             :   ContourFindingBase(ao),
     119           1 :   ones(getPntrToArgument(0)->getRank(),1) {
     120           1 :   if( getPntrToArgument(0)->getRank()<2 ) {
     121           0 :     error("cannot find dividing surface if input grid is one dimensional");
     122             :   }
     123             : 
     124             :   std::string dir;
     125           1 :   parse("SEARCHDIR",dir);
     126           1 :   log.printf("  calculating location of contour on %d dimensional grid \n", getPntrToArgument(0)->getRank()-1 );
     127           1 :   checkRead();
     128             : 
     129             :   Value* gval=getPntrToArgument(0);
     130             :   unsigned n=0;
     131           1 :   gdirs.resize( gval->getRank()-1 );
     132           1 :   gnames.resize( getPntrToArgument(0)->getRank()-1 );
     133             : 
     134           1 :   gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( gval->getPntrToAction() );
     135           1 :   if( !ag ) {
     136           0 :     error("input argument must be a grid");
     137             :   }
     138           2 :   if( getInputGridObject().getGridType()=="fibonacci") {
     139           0 :     error("cannot search for contours in fibonacci grids");
     140             :   }
     141           1 :   std::vector<std::string> argn( ag->getGridCoordinateNames() );
     142             : 
     143           4 :   for(unsigned i=0; i<gval->getRank(); ++i) {
     144           3 :     if( argn[i]==dir ) {
     145           1 :       dir_n=i;
     146             :     } else {
     147           2 :       if( n==gdirs.size() ) {
     148           0 :         error("could not find " + dir + " direction in input grid");
     149             :       }
     150           2 :       gdirs[n]=i;
     151             :       gnames[n]=argn[i];
     152           2 :       n++;
     153             :     }
     154             :   }
     155           1 :   if( n!=(gval->getRank()-1) ) {
     156           0 :     error("output of grid is not understood");
     157             :   }
     158             : 
     159           1 :   std::vector<bool> ipbc( getInputGridObject().getDimension()-1 );
     160           3 :   for(unsigned i=0; i<gdirs.size(); ++i) {
     161           2 :     ipbc[i] = getInputGridObject().isPeriodic(gdirs[i]);
     162             :   }
     163           2 :   gridcoords.setup( "flat", ipbc, 0, 0.0 );
     164             : 
     165             :   // Now add a value
     166           1 :   std::vector<unsigned> shape( getInputGridObject().getDimension()-1 );
     167           1 :   addValueWithDerivatives( shape );
     168           1 :   setNotPeriodic();
     169           1 :   getPntrToComponent(0)->buildDataStore();
     170           2 : }
     171             : 
     172           1 : void FindContourSurface::setupValuesOnFirstStep() {
     173             :   std::vector<double> fspacing;
     174           1 :   std::vector<unsigned> snbins( gridcoords.getDimension() );
     175           1 :   std::vector<std::string> smin( gridcoords.getDimension() ), smax( gridcoords.getDimension() );
     176           3 :   for(unsigned i=0; i<gdirs.size(); ++i) {
     177           4 :     smin[i]=getInputGridObject().getMin()[gdirs[i]];
     178           4 :     smax[i]=getInputGridObject().getMax()[gdirs[i]];
     179           2 :     snbins[i]=getInputGridObject().getNbin(false)[gdirs[i]];
     180             :   }
     181           1 :   gridcoords.setBounds( smin, smax, snbins, fspacing );
     182           2 :   getPntrToComponent(0)->setShape( gridcoords.getNbin(true) );
     183             : 
     184           1 :   std::vector<unsigned> find( gridcoords.getDimension() );
     185           1 :   std::vector<unsigned> ind( gridcoords.getDimension() );
     186         197 :   for(unsigned i=0; i<gridcoords.getNumberOfPoints(); ++i) {
     187         196 :     find.assign( find.size(), 0 );
     188         196 :     gridcoords.getIndices( i, ind );
     189         588 :     for(unsigned j=0; j<gdirs.size(); ++j) {
     190         392 :       find[gdirs[j]]=ind[j];
     191             :     }
     192             :   }
     193             : 
     194             :   // Set the direction in which to look for the contour
     195           1 :   direction.resize( getInputGridObject().getDimension(), 0 );
     196           1 :   direction[dir_n] = 0.999999999*getInputGridObject().getGridSpacing()[dir_n];
     197           2 : }
     198             : 
     199           3 : unsigned FindContourSurface::getNumberOfDerivatives() {
     200           3 :   return gridcoords.getDimension();
     201             : }
     202             : 
     203           3 : std::vector<std::string> FindContourSurface::getGridCoordinateNames() const {
     204           3 :   return gnames;
     205             : }
     206             : 
     207           3 : const gridtools::GridCoordinatesObject& FindContourSurface::getGridCoordinatesObject() const {
     208           3 :   return gridcoords;
     209             : }
     210             : 
     211         588 : void FindContourSurface::performTask( const unsigned& current, MultiValue& myvals ) const {
     212             :   std::vector<unsigned> neighbours;
     213             :   unsigned num_neighbours;
     214             :   unsigned nfound=0;
     215             :   double minv=0, minp;
     216         588 :   std::vector<unsigned> bins_n( getInputGridObject().getNbin(false) );
     217         588 :   unsigned shiftn=current;
     218         588 :   std::vector<unsigned> ind( getInputGridObject().getDimension() );
     219         588 :   std::vector<double> point( getInputGridObject().getDimension() );
     220             : #ifndef DNDEBUG
     221         588 :   std::vector<unsigned> oind( gridcoords.getDimension() );
     222         588 :   gridcoords.getIndices( current, oind );
     223             : #endif
     224       16842 :   for(unsigned i=0; i<bins_n[dir_n]; ++i) {
     225             : #ifndef DNDEBUG
     226       16842 :     std::vector<unsigned> base_ind( getInputGridObject().getDimension() );
     227       16842 :     getInputGridObject().getIndices( shiftn, base_ind );
     228       50526 :     for(unsigned j=0; j<gdirs.size(); ++j) {
     229             :       plumed_dbg_assert( base_ind[gdirs[j]]==oind[j] );
     230             :     }
     231             : #endif
     232             :     // Get the index of the current grid point
     233       16842 :     getInputGridObject().getIndices( shiftn, ind );
     234             :     // Exit if we are at the edge of the grid
     235       16842 :     if( !getInputGridObject().isPeriodic(dir_n) && (ind[dir_n]+1)==bins_n[dir_n] ) {
     236           0 :       shiftn += getInputGridObject().getStride()[dir_n];
     237             :       continue;
     238             :     }
     239             : 
     240             :     // Ensure points with inactive neighbours are ignored
     241       16842 :     getInputGridObject().getNeighbors( ind, ones, num_neighbours, neighbours );
     242             : 
     243             :     // Now get the function value at two points
     244       16842 :     double val1=getPntrToArgument(0)->get( shiftn ) - contour;
     245             :     double val2;
     246       16842 :     if( (ind[dir_n]+1)==bins_n[dir_n] ) {
     247           0 :       val2 = getPntrToArgument(0)->get( current ) - contour;
     248             :     } else {
     249       16842 :       val2=getPntrToArgument(0)->get( shiftn + getInputGridObject().getStride()[dir_n] ) - contour;
     250             :     }
     251             : 
     252             :     // Check if the minimum is bracketed
     253       16842 :     if( val1*val2<0 ) {
     254         588 :       getInputGridObject().getGridPointCoordinates( shiftn, point );
     255         588 :       findContour( direction, point );
     256         588 :       minp=point[dir_n];
     257             :       nfound++;
     258             :       break;
     259             :     }
     260             : 
     261             :     // This moves us on to the next point
     262       16254 :     shiftn += getInputGridObject().getStride()[dir_n];
     263             :   }
     264             :   if( nfound==0 ) {
     265             :     std::string num;
     266           0 :     Tools::convert( getStep(), num );
     267           0 :     error("On step " + num + " failed to find required grid point");
     268             :   }
     269         588 :   myvals.setValue( getConstPntrToComponent(0)->getPositionInStream(), minp );
     270         588 : }
     271             : 
     272         588 : void FindContourSurface::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     273             :     const unsigned& bufstart, std::vector<double>& buffer ) const {
     274             :   plumed_dbg_assert( valindex==0 );
     275         588 :   unsigned istart = bufstart + (1+gridcoords.getDimension())*code;
     276         588 :   unsigned valout = getConstPntrToComponent(0)->getPositionInStream();
     277         588 :   buffer[istart] += myvals.get( valout );
     278         588 : }
     279             : 
     280             : }
     281             : }