Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "FindContour.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "core/PlumedMain.h"
      25             : 
      26             : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR
      27             : /*
      28             : Find an isocontour in a smooth function.
      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 or more dimensions
      34             : it can be 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 values.  In other words, for the function \f$f(x,y)\f$ this action would find a set
      38             : of points \f$\{x_c,y_c\}\f$ that have:
      39             : 
      40             : \f[
      41             : f(x_c,y_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 detected using a variant
      45             : on the marching squares or marching cubes algorithm, which you can find information on here:
      46             : 
      47             : https://en.wikipedia.org/wiki/Marching_squares
      48             : https://en.wikipedia.org/wiki/Marching_cubes
      49             : 
      50             : As such, and unlike \ref FIND_CONTOUR_SURFACE or \ref FIND_SPHERICAL_CONTOUR, the function input to this action can have any dimension.
      51             : Furthermore, the topology of the contour will be determined by the algorithm and does not need to be specified by the user.
      52             : 
      53             : \par Examples
      54             : 
      55             : The input below allows you to calculate something akin to a Willard-Chandler dividing surface \cite wcsurface.
      56             : The simulation cell in this case contains a solid phase and a liquid phase.  The Willard-Chandler surface is the
      57             : surface that separates the parts of the box containing the solid from the parts containing the liquid.  To compute the position
      58             : of this surface  the \ref FCCUBIC symmetry function is calculated for each of the atoms in the system from on the geometry of the
      59             : atoms in the first coordination sphere of each of the atoms.  These quantities are then transformed using a switching function.
      60             : This procedure generates a single number for each atom in the system and this quantity has a value of one for atoms that are in
      61             : parts of the box that resemble the solid structure and zero for atoms that are in parts of the box that resemble the liquid.
      62             : The position of a virtual atom is then computed using \ref CENTER_OF_MULTICOLVAR and a phase field model is constructed using
      63             : \ref MULTICOLVARDENS.  These procedure ensures that we have a continuous function that gives a measure of the average degree of
      64             : solidness at each point in the simulation cell.  The Willard-Chandler dividing surface is calculated by finding a a set of points
      65             : at which the value of this phase field is equal to 0.5.  This set of points is output to file called mycontour.dat.  A new contour
      66             : is found on every single step for each frame that is read in.
      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             : tfcc: MTRANSFORM_MORE DATA=fcc LOWMEM SWITCH={SMAP R_0=0.5 A=8 B=8}
      76             : center: CENTER_OF_MULTICOLVAR DATA=tfcc
      77             : 
      78             : dens: MULTICOLVARDENS ...
      79             :   DATA=tfcc ORIGIN=center DIR=xyz
      80             :   NBINS=80,80,80 BANDWIDTH=1.0,1.0,1.0 STRIDE=1 CLEAR=1
      81             : ...
      82             : 
      83             : FIND_CONTOUR GRID=dens CONTOUR=0.5 FILE=mycontour.xyz
      84             : \endplumedfile
      85             : 
      86             : */
      87             : //+ENDPLUMEDOC
      88             : 
      89             : namespace PLMD {
      90             : namespace contour {
      91             : 
      92             : PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
      93             : 
      94           5 : void FindContour::registerKeywords( Keywords& keys ) {
      95           5 :   ContourFindingBase::registerKeywords( keys ); ActionWithValue::useCustomisableComponents(keys);
      96             : // We want a better way of doing this bit
      97          10 :   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");
      98           5 : }
      99             : 
     100           2 : FindContour::FindContour(const ActionOptions&ao):
     101             :   Action(ao),
     102           2 :   ContourFindingBase(ao)
     103             : {
     104           2 :   parse("BUFFER",gbuffer);
     105           2 :   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);
     106           2 :   checkRead();
     107             : 
     108           2 :   gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
     109           2 :   std::vector<std::string> argn( ag->getGridCoordinateNames() );
     110             : 
     111           2 :   std::vector<unsigned> shape(1); shape[0]=0;
     112           8 :   for(unsigned i=0; i<argn.size(); ++i ) {
     113           6 :     addComponent( argn[i], shape ); componentIsNotPeriodic( argn[i] ); getPntrToComponent(i)->buildDataStore();
     114             :   }
     115             :   // Check for task reduction
     116           2 :   updateTaskListReductionStatus();
     117           2 : }
     118             : 
     119           3 : std::string FindContour::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
     120           6 :   return "a vector of coordinates for the contour along the " + cname + " direction";
     121             : }
     122             : 
     123           1 : void FindContour::setupValuesOnFirstStep() {
     124           1 :   std::vector<unsigned> shape(1); shape[0] = getPntrToArgument(0)->getRank()*getPntrToArgument(0)->getNumberOfValues();
     125           4 :   for(unsigned i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->setShape( shape );
     126           1 :   active_cells.resize( shape[0] );
     127           1 : }
     128             : 
     129           0 : unsigned FindContour::getNumberOfDerivatives() {
     130           0 :   return 0;
     131             : }
     132             : 
     133           2 : void FindContour::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
     134           2 :   task_reducing_actions.push_back(this);
     135           2 : }
     136             : 
     137           2 : void FindContour::getNumberOfTasks( unsigned& ntasks ) {
     138           2 :   ntasks = active_cells.size();
     139             : 
     140             :   Value* gval=getPntrToArgument(0);
     141           2 :   unsigned npoints = gval->getNumberOfValues();
     142           2 :   std::vector<unsigned> ind( gval->getRank() );
     143           2 :   std::vector<unsigned> ones( gval->getRank(), 1 );
     144           2 :   std::vector<unsigned> nbin( getInputGridObject().getNbin( false ) );
     145             :   unsigned num_neighbours; std::vector<unsigned> neighbours;
     146             : 
     147             :   std::fill( active_cells.begin(), active_cells.end(), 0 );
     148       10978 :   for(unsigned i=0; i<npoints; ++i) {
     149             :     // Get the index of the current grid point
     150       10976 :     getInputGridObject().getIndices( i, ind );
     151       10976 :     getInputGridObject().getNeighbors( ind, ones, num_neighbours, neighbours );
     152             :     // Get the value of a point on the grid
     153       10976 :     double val1=gval->get( i ) - contour;
     154             :     bool edge=false;
     155       43904 :     for(unsigned j=0; j<gval->getRank(); ++j) {
     156             :       // Make sure we don't search at the edge of the grid
     157       32928 :       if( !getInputGridObject().isPeriodic(j) && (ind[j]+1)==nbin[j] ) continue;
     158       32928 :       else if( (ind[j]+1)==nbin[j] ) { edge=true; ind[j]=0; }
     159       30968 :       else ind[j]+=1;
     160       32928 :       double val2=gval->get( getInputGridObject().getIndex(ind) ) - contour;
     161       32928 :       if( val1*val2<0 ) active_cells[gval->getRank()*i + j] = 1;
     162       32928 :       if( getInputGridObject().isPeriodic(j) && edge ) { edge=false; ind[j]=nbin[j]-1; }
     163       30968 :       else ind[j]-=1;
     164             :     }
     165             :   }
     166           2 : }
     167             : 
     168       32928 : int FindContour::checkTaskStatus( const unsigned& taskno, int& flag ) const {
     169       32928 :   if( active_cells[taskno]>0 ) return 1;
     170             :   return 0;
     171             : }
     172             : 
     173        1108 : void FindContour::performTask( const unsigned& current, MultiValue& myvals ) const {
     174             :   // Retrieve the initial grid point coordinates
     175        1108 :   unsigned gpoint = std::floor( current / getPntrToArgument(0)->getRank() );
     176        1108 :   std::vector<double> point( getPntrToArgument(0)->getRank() );
     177        1108 :   getInputGridObject().getGridPointCoordinates( gpoint, point );
     178             : 
     179             :   // Retrieve the direction we are searching for the contour
     180        1108 :   unsigned gdir = current%(getPntrToArgument(0)->getRank() );
     181        1108 :   std::vector<double> direction( getPntrToArgument(0)->getRank(), 0 );
     182        1108 :   direction[gdir] = 0.999999999*getInputGridObject().getGridSpacing()[gdir];
     183             : 
     184             :   // Now find the contour
     185             :   findContour( direction, point );
     186             :   // And transfer to the store data vessel
     187        4432 :   for(unsigned i=0; i<getPntrToArgument(0)->getRank(); ++i) myvals.setValue( getConstPntrToComponent(i)->getPositionInStream(), point[i] );
     188        1108 : }
     189             : 
     190             : }
     191             : }