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             : 
      25             : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR_SURFACE
      26             : /*
      27             : Find an isocontour by searching along either the x, y or z direction.
      28             : 
      29             : As discussed in the part of the manual on \ref Analysis PLUMED contains a number of tools that allow you to calculate
      30             : a function on a grid.  The function on this grid might be a \ref HISTOGRAM as a function of a few collective variables
      31             : or it might be a phase field that has been calculated using \ref MULTICOLVARDENS.  If this function has one or two input
      32             : arguments it is relatively straightforward to plot the function.  If by contrast the data has a three dimensions it can be
      33             : difficult to visualize.
      34             : 
      35             : This action provides one tool for visualizing these functions.  It can be used to search for a set of points on a contour
      36             : 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
      37             : of points \f$\{x_c,y_c,z_c\}\f$ that have:
      38             : 
      39             : \f[
      40             : f(x_c,y_c,z_c) - c = 0
      41             : \f]
      42             : 
      43             : 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
      44             : 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
      45             : function evaluated on a grid that gives us the height of the interface as a function of two coordinates.
      46             : 
      47             : 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
      48             : 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
      49             : function have the appropriate topology you should use \ref FIND_CONTOUR in place of \ref FIND_CONTOUR_SURFACE.
      50             : 
      51             : 
      52             : \par Examples
      53             : 
      54             : 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
      55             : 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
      56             : 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
      57             : 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
      58             : 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
      59             : \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
      60             : 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
      61             : 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
      62             : 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$
      63             : position.  This grid is then output to a file called contour2.dat.
      64             : 
      65             : 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
      66             : on every step.
      67             : 
      68             : \plumedfile
      69             : UNITS NATURAL
      70             : FCCUBIC ...
      71             :   SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
      72             :   ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619 LABEL=fcc
      73             : ... FCCUBIC
      74             : 
      75             : 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
      76             : 
      77             : ss2: FIND_CONTOUR_SURFACE GRID=dens2 CONTOUR=0.42 SEARCHDIR=z STRIDE=1 CLEAR=1
      78             : DUMPGRID GRID=ss2 FILE=contour2.dat FMT=%8.4f STRIDE=1
      79             : \endplumedfile
      80             : 
      81             : */
      82             : //+ENDPLUMEDOC
      83             : 
      84             : namespace PLMD {
      85             : namespace gridtools {
      86             : 
      87             : class FindContourSurface : public ContourFindingBase {
      88             : private:
      89             :   bool firsttime;
      90             :   unsigned dir_n;
      91             :   unsigned gbuffer;
      92             :   std::vector<unsigned> ones;
      93             :   std::vector<unsigned> gdirs;
      94             :   std::vector<double> direction;
      95             : public:
      96             :   static void registerKeywords( Keywords& keys );
      97             :   explicit FindContourSurface(const ActionOptions&ao);
      98          12 :   unsigned getNumberOfQuantities() const override { return 2; }
      99           0 :   bool checkAllActive() const override { return gbuffer==0; }
     100             :   void clearAverage() override;
     101             :   void prepareForAveraging() override;
     102             :   void compute( const unsigned& current, MultiValue& myvals ) const override;
     103             :   void finishAveraging() override;
     104             : };
     105             : 
     106       10421 : PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
     107             : 
     108           2 : void FindContourSurface::registerKeywords( Keywords& keys ) {
     109           2 :   ContourFindingBase::registerKeywords( keys );
     110           4 :   keys.add("compulsory","SEARCHDIR","In which directions do you wish to search for the contour.");
     111           4 :   keys.add("compulsory","BUFFER","0","number of buffer grid points around location where grid was found on last step.  If this is zero the full grid is calculated on each step");
     112           2 : }
     113             : 
     114           1 : FindContourSurface::FindContourSurface(const ActionOptions&ao):
     115             :   Action(ao),
     116             :   ContourFindingBase(ao),
     117           1 :   firsttime(true),
     118           1 :   ones(ingrid->getDimension(),1)
     119             : {
     120           1 :   if( ingrid->getDimension()<2 ) error("cannot find dividing surface if input grid is one dimensional");
     121             : 
     122           2 :   std::string dir; parse("SEARCHDIR",dir); parse("BUFFER",gbuffer);
     123           1 :   log.printf("  calculating location of contour on %d dimensional grid \n", ingrid->getDimension()-1 );
     124           1 :   if( gbuffer>0 ) log.printf("  after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
     125           1 :   checkRead();
     126             : 
     127           1 :   unsigned n=0; gdirs.resize( ingrid->getDimension()-1 );
     128           4 :   for(unsigned i=0; i<ingrid->getDimension(); ++i) {
     129           3 :     if( ingrid->getComponentName(i)==dir ) {
     130           1 :       dir_n=i;
     131             :     } else {
     132           2 :       if( n==gdirs.size() ) error("could not find " + dir + " direction in input grid");
     133           2 :       gdirs[n]=i; n++;
     134             :     }
     135             :   }
     136           1 :   if( n!=(ingrid->getDimension()-1) ) error("output of grid is not understood");
     137             : 
     138             :   // Create the input from the old string
     139           2 :   std::string vstring = "COMPONENTS=" + getLabel() + " COORDINATES=" + ingrid->getComponentName( gdirs[0] );
     140           2 :   for(unsigned i=1; i<gdirs.size(); ++i) vstring += "," + ingrid->getComponentName( gdirs[i] );
     141             :   vstring += " PBC=";
     142           1 :   if( ingrid->isPeriodic(gdirs[0]) ) vstring+="T";
     143             :   else vstring+="F";
     144           2 :   for(unsigned i=1; i<gdirs.size(); ++i) {
     145           1 :     if( ingrid->isPeriodic(gdirs[i]) ) vstring+=",T"; else vstring+=",F";
     146             :   }
     147           2 :   auto grid=createGrid( "grid", vstring ); grid->setNoDerivatives();
     148           1 :   setAveragingAction( std::move(grid), true );
     149           2 : }
     150             : 
     151           3 : void FindContourSurface::clearAverage() {
     152             :   // Set the boundaries of the output grid
     153           3 :   std::vector<double> fspacing; std::vector<unsigned> snbins( ingrid->getDimension()-1 );
     154           3 :   std::vector<std::string> smin( ingrid->getDimension()-1 ), smax( ingrid->getDimension()-1 );
     155           9 :   for(unsigned i=0; i<gdirs.size(); ++i) {
     156          18 :     smin[i]=ingrid->getMin()[gdirs[i]]; smax[i]=ingrid->getMax()[gdirs[i]];
     157          12 :     snbins[i]=ingrid->getNbin()[gdirs[i]];
     158             :   }
     159           3 :   mygrid->setBounds( smin, smax, snbins, fspacing); resizeFunctions();
     160           3 :   ActionWithAveraging::clearAverage();
     161           6 : }
     162             : 
     163           3 : void FindContourSurface::prepareForAveraging() {
     164             :   // Create a task list if first time
     165           3 :   if( firsttime ) {
     166           1 :     std::vector<unsigned> find( ingrid->getDimension() );
     167           1 :     std::vector<unsigned> ind( mygrid->getDimension() );
     168         197 :     for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
     169         196 :       find.assign( find.size(), 0 ); mygrid->getIndices( i, ind );
     170         588 :       for(unsigned j=0; j<gdirs.size(); ++j) find[gdirs[j]]=ind[j];
     171             :       // Current will be set equal to the start point for this grid index
     172         196 :       addTaskToList( ingrid->getIndex(find) );
     173             :     }
     174             :     // And prepare the task list
     175           1 :     deactivateAllTasks();
     176         197 :     for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
     177           1 :     lockContributors();
     178             :     // Set the direction in which to look for the contour
     179           1 :     direction.resize( ingrid->getDimension(), 0 );
     180           1 :     direction[dir_n] = 0.999999999*ingrid->getGridSpacing()[dir_n];
     181             :   }
     182           3 : }
     183             : 
     184           3 : void FindContourSurface::finishAveraging() {
     185             :   ContourFindingBase::finishAveraging();
     186             :   // And update the list of active grid points
     187           3 :   if( gbuffer>0 ) {
     188           3 :     std::vector<double> dx( ingrid->getGridSpacing() );
     189           3 :     std::vector<double> point( ingrid->getDimension() );
     190           3 :     std::vector<double> lpoint( mygrid->getDimension() );
     191             :     std::vector<unsigned> neighbours; unsigned num_neighbours;
     192           3 :     std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
     193           3 :     std::vector<bool> active( ingrid->getNumberOfPoints(), false );
     194           3 :     std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
     195         591 :     for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
     196             :       // Retrieve the coordinates of this grid point
     197         588 :       mygrid->getGridPointCoordinates( i, lpoint );
     198         588 :       point[dir_n] = mygrid->getGridElement( i, 0 );
     199             :       // 0.5*dx added here to prevent problems with flooring of grid points
     200        1764 :       for(unsigned j=0; j<gdirs.size(); ++j) point[gdirs[j]]=lpoint[j] + 0.5*dx[gdirs[j]];
     201         588 :       ingrid->getIndices( point, ugrid_indices );
     202             :       // Now activate buffer region
     203         588 :       ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
     204       16464 :       for(unsigned n=0; n<num_neighbours; ++n) active[ neighbours[n] ]=true;
     205             :     }
     206           3 :     ingrid->activateThesePoints( active );
     207             :   }
     208           3 :   firsttime=false;
     209           3 : }
     210             : 
     211         588 : void FindContourSurface::compute( const unsigned& current, MultiValue& myvals ) const {
     212             :   std::vector<unsigned> neighbours; unsigned num_neighbours; unsigned nfound=0; double minp=0;
     213         588 :   std::vector<unsigned> bins_n( ingrid->getNbin() ); unsigned shiftn=current;
     214         588 :   std::vector<unsigned> ind( ingrid->getDimension() ); std::vector<double> point( ingrid->getDimension() );
     215             : #ifndef NDEBUG
     216             :   std::vector<unsigned> oind( mygrid->getDimension() ); mygrid->getIndices( current, oind );
     217             : #endif
     218       16749 :   for(unsigned i=0; i<bins_n[dir_n]; ++i) {
     219             : #ifndef NDEBUG
     220             :     std::vector<unsigned> base_ind( ingrid->getDimension() ); ingrid->getIndices( shiftn, base_ind );
     221             :     for(unsigned j=0; j<gdirs.size(); ++j) plumed_dbg_assert( base_ind[gdirs[j]]==oind[j] );
     222             : #endif
     223             :     // Ensure inactive grid points are ignored
     224       16749 :     if( ingrid->inactive( shiftn ) ) { shiftn += ingrid->getStride()[dir_n]; continue; }
     225             :     // Get the index of the current grid point
     226        7857 :     ingrid->getIndices( shiftn, ind );
     227             :     // Exit if we are at the edge of the grid
     228        7857 :     if( !ingrid->isPeriodic(dir_n) && (ind[dir_n]+1)==bins_n[dir_n] ) {
     229           0 :       shiftn += ingrid->getStride()[dir_n]; continue;
     230             :     }
     231             : 
     232             :     // Ensure points with inactive neighbours are ignored
     233        7857 :     ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
     234             :     bool cycle=false;
     235      179020 :     for(unsigned j=0; j<num_neighbours; ++j) {
     236      172753 :       if( ingrid->inactive( neighbours[j]) ) { cycle=true; break; }
     237             :     }
     238        7857 :     if( cycle ) { shiftn += ingrid->getStride()[dir_n]; continue; }
     239             : 
     240             :     // Now get the function value at two points
     241        6267 :     double val1=getFunctionValue( shiftn ) - contour; double val2;
     242        6267 :     if( (ind[dir_n]+1)==bins_n[dir_n] ) val2 = getFunctionValue( current ) - contour;
     243        6267 :     else val2=getFunctionValue( shiftn + ingrid->getStride()[dir_n] ) - contour;
     244             : 
     245             :     // Check if the minimum is bracketed
     246        6267 :     if( val1*val2<0 ) {
     247         588 :       ingrid->getGridPointCoordinates( shiftn, point ); findContour( direction, point );
     248         588 :       minp=point[dir_n]; nfound++; break;
     249             :     }
     250             : 
     251             : 
     252             :     // This moves us on to the next point
     253        5679 :     shiftn += ingrid->getStride()[dir_n];
     254             :   }
     255             :   if( nfound==0 ) {
     256           0 :     std::string num; Tools::convert( getStep(), num );
     257           0 :     error("On step " + num + " failed to find required grid point");
     258             :   }
     259             :   myvals.setValue( 1, minp );
     260         588 : }
     261             : 
     262             : }
     263             : }