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 );
      96           5 :   ActionWithValue::useCustomisableComponents(keys);
      97             : // We want a better way of doing this bit
      98           5 :   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");
      99           5 : }
     100             : 
     101           2 : FindContour::FindContour(const ActionOptions&ao):
     102             :   Action(ao),
     103           2 :   ContourFindingBase(ao) {
     104           2 :   parse("BUFFER",gbuffer);
     105           2 :   if( gbuffer>0 ) {
     106           0 :     log.printf("  after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
     107             :   }
     108           2 :   checkRead();
     109             : 
     110           2 :   gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
     111           2 :   std::vector<std::string> argn( ag->getGridCoordinateNames() );
     112             : 
     113           2 :   std::vector<unsigned> shape(1);
     114           2 :   shape[0]=0;
     115           8 :   for(unsigned i=0; i<argn.size(); ++i ) {
     116           6 :     addComponent( argn[i], shape );
     117           6 :     componentIsNotPeriodic( argn[i] );
     118           6 :     getPntrToComponent(i)->buildDataStore();
     119             :   }
     120             :   // Check for task reduction
     121           2 :   updateTaskListReductionStatus();
     122           2 : }
     123             : 
     124           3 : std::string FindContour::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
     125           6 :   return "a vector of coordinates for the contour along the " + cname + " direction";
     126             : }
     127             : 
     128           1 : void FindContour::setupValuesOnFirstStep() {
     129           1 :   std::vector<unsigned> shape(1);
     130           1 :   shape[0] = getPntrToArgument(0)->getRank()*getPntrToArgument(0)->getNumberOfValues();
     131           4 :   for(unsigned i=0; i<getNumberOfComponents(); ++i) {
     132           3 :     getPntrToComponent(i)->setShape( shape );
     133             :   }
     134           1 :   active_cells.resize( shape[0] );
     135           1 : }
     136             : 
     137           0 : unsigned FindContour::getNumberOfDerivatives() {
     138           0 :   return 0;
     139             : }
     140             : 
     141           2 : void FindContour::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
     142           2 :   task_reducing_actions.push_back(this);
     143           2 : }
     144             : 
     145           2 : void FindContour::getNumberOfTasks( unsigned& ntasks ) {
     146           2 :   ntasks = active_cells.size();
     147             : 
     148             :   Value* gval=getPntrToArgument(0);
     149           2 :   unsigned npoints = gval->getNumberOfValues();
     150           2 :   std::vector<unsigned> ind( gval->getRank() );
     151           2 :   std::vector<unsigned> ones( gval->getRank(), 1 );
     152           2 :   std::vector<unsigned> nbin( getInputGridObject().getNbin( false ) );
     153             :   unsigned num_neighbours;
     154             :   std::vector<unsigned> neighbours;
     155             : 
     156             :   std::fill( active_cells.begin(), active_cells.end(), 0 );
     157       10978 :   for(unsigned i=0; i<npoints; ++i) {
     158             :     // Get the index of the current grid point
     159       10976 :     getInputGridObject().getIndices( i, ind );
     160       10976 :     getInputGridObject().getNeighbors( ind, ones, num_neighbours, neighbours );
     161             :     // Get the value of a point on the grid
     162       10976 :     double val1=gval->get( i ) - contour;
     163             :     bool edge=false;
     164       43904 :     for(unsigned j=0; j<gval->getRank(); ++j) {
     165             :       // Make sure we don't search at the edge of the grid
     166       32928 :       if( !getInputGridObject().isPeriodic(j) && (ind[j]+1)==nbin[j] ) {
     167           0 :         continue;
     168       32928 :       } else if( (ind[j]+1)==nbin[j] ) {
     169             :         edge=true;
     170        1960 :         ind[j]=0;
     171             :       } else {
     172       30968 :         ind[j]+=1;
     173             :       }
     174       32928 :       double val2=gval->get( getInputGridObject().getIndex(ind) ) - contour;
     175       32928 :       if( val1*val2<0 ) {
     176        1108 :         active_cells[gval->getRank()*i + j] = 1;
     177             :       }
     178       32928 :       if( getInputGridObject().isPeriodic(j) && edge ) {
     179             :         edge=false;
     180        1960 :         ind[j]=nbin[j]-1;
     181             :       } else {
     182       30968 :         ind[j]-=1;
     183             :       }
     184             :     }
     185             :   }
     186           2 : }
     187             : 
     188       32928 : int FindContour::checkTaskStatus( const unsigned& taskno, int& flag ) const {
     189       32928 :   if( active_cells[taskno]>0 ) {
     190        1108 :     return 1;
     191             :   }
     192             :   return 0;
     193             : }
     194             : 
     195        1108 : void FindContour::performTask( const unsigned& current, MultiValue& myvals ) const {
     196             :   // Retrieve the initial grid point coordinates
     197        1108 :   unsigned gpoint = std::floor( current / getPntrToArgument(0)->getRank() );
     198        1108 :   std::vector<double> point( getPntrToArgument(0)->getRank() );
     199        1108 :   getInputGridObject().getGridPointCoordinates( gpoint, point );
     200             : 
     201             :   // Retrieve the direction we are searching for the contour
     202        1108 :   unsigned gdir = current%(getPntrToArgument(0)->getRank() );
     203        1108 :   std::vector<double> direction( getPntrToArgument(0)->getRank(), 0 );
     204        1108 :   direction[gdir] = 0.999999999*getInputGridObject().getGridSpacing()[gdir];
     205             : 
     206             :   // Now find the contour
     207             :   findContour( direction, point );
     208             :   // And transfer to the store data vessel
     209        4432 :   for(unsigned i=0; i<getPntrToArgument(0)->getRank(); ++i) {
     210        3324 :     myvals.setValue( getConstPntrToComponent(i)->getPositionInStream(), point[i] );
     211             :   }
     212        1108 : }
     213             : 
     214             : }
     215             : }