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           6 :   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             : {
     121           1 :   if( getPntrToArgument(0)->getRank()<2 ) error("cannot find dividing surface if input grid is one dimensional");
     122             : 
     123           1 :   std::string dir; parse("SEARCHDIR",dir);
     124           1 :   log.printf("  calculating location of contour on %d dimensional grid \n", getPntrToArgument(0)->getRank()-1 );
     125           1 :   checkRead();
     126             : 
     127             :   Value* gval=getPntrToArgument(0); unsigned n=0;
     128           1 :   gdirs.resize( gval->getRank()-1 ); gnames.resize( getPntrToArgument(0)->getRank()-1 );
     129             : 
     130           1 :   gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( gval->getPntrToAction() );
     131           1 :   if( !ag ) error("input argument must be a grid");
     132           2 :   if( getInputGridObject().getGridType()=="fibonacci") error("cannot search for contours in fibonacci grids");
     133           1 :   std::vector<std::string> argn( ag->getGridCoordinateNames() );
     134             : 
     135           4 :   for(unsigned i=0; i<gval->getRank(); ++i) {
     136           3 :     if( argn[i]==dir ) {
     137           1 :       dir_n=i;
     138             :     } else {
     139           2 :       if( n==gdirs.size() ) error("could not find " + dir + " direction in input grid");
     140           4 :       gdirs[n]=i; gnames[n]=argn[i]; n++;
     141             :     }
     142             :   }
     143           1 :   if( n!=(gval->getRank()-1) ) error("output of grid is not understood");
     144             : 
     145           1 :   std::vector<bool> ipbc( getInputGridObject().getDimension()-1 );
     146           3 :   for(unsigned i=0; i<gdirs.size(); ++i) ipbc[i] = getInputGridObject().isPeriodic(gdirs[i]);
     147           2 :   gridcoords.setup( "flat", ipbc, 0, 0.0 );
     148             : 
     149             :   // Now add a value
     150           1 :   std::vector<unsigned> shape( getInputGridObject().getDimension()-1 );
     151           1 :   addValueWithDerivatives( shape ); setNotPeriodic();
     152           1 :   getPntrToComponent(0)->buildDataStore();
     153           2 : }
     154             : 
     155           1 : void FindContourSurface::setupValuesOnFirstStep() {
     156           1 :   std::vector<double> fspacing; std::vector<unsigned> snbins( gridcoords.getDimension() );
     157           1 :   std::vector<std::string> smin( gridcoords.getDimension() ), smax( gridcoords.getDimension() );
     158           3 :   for(unsigned i=0; i<gdirs.size(); ++i) {
     159           6 :     smin[i]=getInputGridObject().getMin()[gdirs[i]]; smax[i]=getInputGridObject().getMax()[gdirs[i]];
     160           2 :     snbins[i]=getInputGridObject().getNbin(false)[gdirs[i]];
     161             :   }
     162           1 :   gridcoords.setBounds( smin, smax, snbins, fspacing );
     163           2 :   getPntrToComponent(0)->setShape( gridcoords.getNbin(true) );
     164             : 
     165           1 :   std::vector<unsigned> find( gridcoords.getDimension() );
     166           1 :   std::vector<unsigned> ind( gridcoords.getDimension() );
     167         197 :   for(unsigned i=0; i<gridcoords.getNumberOfPoints(); ++i) {
     168         196 :     find.assign( find.size(), 0 ); gridcoords.getIndices( i, ind );
     169         588 :     for(unsigned j=0; j<gdirs.size(); ++j) find[gdirs[j]]=ind[j];
     170             :   }
     171             : 
     172             :   // Set the direction in which to look for the contour
     173           1 :   direction.resize( getInputGridObject().getDimension(), 0 );
     174           1 :   direction[dir_n] = 0.999999999*getInputGridObject().getGridSpacing()[dir_n];
     175           2 : }
     176             : 
     177           3 : unsigned FindContourSurface::getNumberOfDerivatives() {
     178           3 :   return gridcoords.getDimension();
     179             : }
     180             : 
     181           3 : std::vector<std::string> FindContourSurface::getGridCoordinateNames() const {
     182           3 :   return gnames;
     183             : }
     184             : 
     185           3 : const gridtools::GridCoordinatesObject& FindContourSurface::getGridCoordinatesObject() const {
     186           3 :   return gridcoords;
     187             : }
     188             : 
     189         588 : void FindContourSurface::performTask( const unsigned& current, MultiValue& myvals ) const {
     190             :   std::vector<unsigned> neighbours; unsigned num_neighbours; unsigned nfound=0; double minv=0, minp;
     191         588 :   std::vector<unsigned> bins_n( getInputGridObject().getNbin(false) ); unsigned shiftn=current;
     192         588 :   std::vector<unsigned> ind( getInputGridObject().getDimension() ); std::vector<double> point( getInputGridObject().getDimension() );
     193             : #ifndef DNDEBUG
     194         588 :   std::vector<unsigned> oind( gridcoords.getDimension() ); gridcoords.getIndices( current, oind );
     195             : #endif
     196       16842 :   for(unsigned i=0; i<bins_n[dir_n]; ++i) {
     197             : #ifndef DNDEBUG
     198       16842 :     std::vector<unsigned> base_ind( getInputGridObject().getDimension() ); getInputGridObject().getIndices( shiftn, base_ind );
     199       50526 :     for(unsigned j=0; j<gdirs.size(); ++j) plumed_dbg_assert( base_ind[gdirs[j]]==oind[j] );
     200             : #endif
     201             :     // Get the index of the current grid point
     202       16842 :     getInputGridObject().getIndices( shiftn, ind );
     203             :     // Exit if we are at the edge of the grid
     204       16842 :     if( !getInputGridObject().isPeriodic(dir_n) && (ind[dir_n]+1)==bins_n[dir_n] ) {
     205           0 :       shiftn += getInputGridObject().getStride()[dir_n]; continue;
     206             :     }
     207             : 
     208             :     // Ensure points with inactive neighbours are ignored
     209       16842 :     getInputGridObject().getNeighbors( ind, ones, num_neighbours, neighbours );
     210             : 
     211             :     // Now get the function value at two points
     212       16842 :     double val1=getPntrToArgument(0)->get( shiftn ) - contour; double val2;
     213       16842 :     if( (ind[dir_n]+1)==bins_n[dir_n] ) val2 = getPntrToArgument(0)->get( current ) - contour;
     214       16842 :     else val2=getPntrToArgument(0)->get( shiftn + getInputGridObject().getStride()[dir_n] ) - contour;
     215             : 
     216             :     // Check if the minimum is bracketed
     217       16842 :     if( val1*val2<0 ) {
     218         588 :       getInputGridObject().getGridPointCoordinates( shiftn, point ); findContour( direction, point );
     219         588 :       minp=point[dir_n]; nfound++; break;
     220             :     }
     221             : 
     222             :     // This moves us on to the next point
     223       16254 :     shiftn += getInputGridObject().getStride()[dir_n];
     224             :   }
     225             :   if( nfound==0 ) {
     226           0 :     std::string num; Tools::convert( getStep(), num );
     227           0 :     error("On step " + num + " failed to find required grid point");
     228             :   }
     229         588 :   myvals.setValue( getConstPntrToComponent(0)->getPositionInStream(), minp );
     230         588 : }
     231             : 
     232         588 : void FindContourSurface::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     233             :     const unsigned& bufstart, std::vector<double>& buffer ) const {
     234         588 :   plumed_dbg_assert( valindex==0 ); unsigned istart = bufstart + (1+gridcoords.getDimension())*code;
     235         588 :   unsigned valout = getConstPntrToComponent(0)->getPositionInStream(); buffer[istart] += myvals.get( valout );
     236         588 : }
     237             : 
     238             : }
     239             : }